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ABSTRACT 

Context. Possible effects of magnetic fields in core collapse supernovae rely on an efficient amplification of the weak pre-collapse 
fields. It has been suggested that the magneto-rotational instability (MRl) leads to a rapid growth for these weak seed fields. Although 
plenty of MRl studies exist for accretion disks, the application of their results to core collapse supernovae is inhibited as the physics 
of supernova cores is substantially different from that of accretion discs. 

Aims. We address the problem of growth and saturation of the MRl in core collapse supernovae by studying its evolution by means 
of semi-global simulations, which combine elements of global and local simulations by taking the presence of global background 
gradients into account and using a local computational grid. We investigate, in particular, the termination of the growth of the MRl 
and the properties of the turbulence in the saturated state. 

Methods. We analyze the dispersion relation of the MRl to identify different regimes of the instability. This analysis is complemented 
by semi-global ideal MHD simulations, where we consider core matter in a local computational box (size ~ 1 km) rotating at sub- 
Keplerian velocity and where we allow for the presence of a radial entropy gradient, but neglect neutrino radiation. 
Results. We identify six regimes of the MRl depending on the ratio of the entropy and angular velocity gradient. Our numerical models 
confirm the instability criteria and growth rates for all regimes relevant to core-collapse supernovae. The MRl grows exponentially on 
time scales of milliseconds, the flow and magnetic field geometries being dominated by channel flows. We find MHD turbulence and 
efficient transport of angular momentum. The MRl growth ceases once the channels are disrupted by resistive instabilities (stemming 
from to the finite conductivity of the numerical code), and MHD turbulence sets in. From an analysis of the growth rates of the 
resistive instabilities, we deduce scaling laws for the termination amplitude of the MRl, which agree well with our numerical models. 
We determine the dependence of the development of large-scale coherent flow structures in the saturated state on the aspect ratio of 
the simulation boxes. 

Conclusions. The MRl can grow rapidly under the conditions considered here, i.e., a rapidly rotating core in hydrostatic equilibrium, 
possibly endowed with a nonvanishing entropy gradient, leading to fields exceeding 10" G. More investigations are required to cover 
the parameter space more comprehensively and to include more physical effects. 
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1. Introduction mentum redistribution and for the conversion of rotational en- 
ergy into thermal energy of the gas. Imparting additional ther- 

The magneto-rotational instability (MRl) ('Balbu s & Hawleyl mal energy into the post-shock stellar matter the MRl might 

il991) is a local linear instability of weakly magnetized differ- thus be important for the currendy favored neu trino-d rive n core 

entially rotating fluids. A large number of analytic, as well as collapse supernova expl osion mechanism (e.g. iThompson et al.l 

numerical studies support the assertion that the MRl is the main l2005t iJanka et al] l2007i) . although possibly only for rapidly 

agent for exciting turbulence in Keplerian accretion disks (for and strongly differentially rotating progenitors. Furthermore, 

a review, see, e.g.. lBalbus & Hawlevlll998 ). The MRl amplifies the growth of the magnetic field resulting from the MRl may 

seed perturbations exponentially with time until turbulence sets provide the adequate physical conditions in the collapsed core 

in. In the turbulent state, the magnetic field, b, gives rise to a to launch bipolar outflows, w hich result in gamma-ray bursts 

non-vanishing (spatial and temporal) mean Maxwell stress ten- (lAlov & Obergaulingej2007 ). As the physical conditions in ac- 

sor Mij - bibj. Simulations of accretion disks show a high neg- cretion disks and stars differ significandy, and as only a few an- 

ative mean value of the component M^j^ (where vj and <p are the alytic studies of the MRl in stars exist (e.g., IAchesonlll978h . it 

radial and azimuthal coordinate of a cylindrical coordinate sys- remains unclear whether existing results on the MRl in disks ap- 

tem), which gives rise to an efficient outward transport of angular ply to stai-s, and particularly to supernovae, as well. 

mo mentum. Numerical simulations of the MRl face a severe problem, 

lAkivama et al.l (l2003h pointed out that the layers surround- since the growth rate of MRI-unstable modes depends on the 

ing the nascent proto-neutron star quite generically fulfill the product of the initial field strength and the wave number of the 

MRl instability criteria. Consequently, any (weak) seed mag- mode. For a weak field, only fairly short modes grow rapidly, 

netic field will be amplified exponentially. In the saturated state Simulations of astrophysical flows, on the other hand, often fail 

of the MRl instability, sustained magneto-hydrodynamic turbu- to resolve just those modes, as it would require prohibitively 

lence might then provide an efficient means for an angular mo- high computational costs to cover spatial scales ranging from 
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the global extent of the astrophysical system (which may be 
much larger than the MRI-unstable region) down to the wave- 
lengths of the fastest growing MRI modes. This dependence 
of the growth rate on the wavelength of the mode suggests a 
twofold approximate numerical approach: one either performs 
simulations which properly cover the global scales of the astro- 
physical system foregoing to resolve the small scales set by the 
wavelengths of the fastest growing MRI modes (global simula- 
tions), or vice versa (local simulations). 

Local simulations evolve only a small part of the entire MRI- 
unstable system, known as the shearing box. However, infor- 
mation on the scales exceeding the size of the computational 
grid has to be provided using suitably chosen boundary condi- 
tions. No unique recipe exists for this procedure, but the use 
of reflecting and periodic boundaries is a common practice. 
In most studies of accretion disks, the boundary conditions of 
lHawlev&Ba"ibusldl992h are used in radial direction, which are 
essentially periodic boundary conditions but account also for the 
relative shear between the inner and the outer radial edge of the 
grid. They are often combined with a Galilei transformation into 
a frame of reference co-rotating at the mean angular velocity 
of the shearing box, and a linearizion of the angular velocity 
within the box. Local (shearing box) simulations using this kind 
of boundary treatment are commonly called shearing-sheet sim- 
ulations. 

The general drawback of local simulations obviously lies in 
their inability to account accurately for large-scale phenomena. 
In addition, there is only a limited possibility to model global 
gradients other than differential rotation in shearing-sheet sim- 
ulations. Independent of the boundary treatment only modes 
with a wavelength less than the size of the grid can be ex- 
cited, i.e. modes with a wavelength comparable to the dimen- 
sions of the whole system cannot develop. Consequently, MRI- 
driven turbulence may saturate at a level determined (at least par- 
tially) by numerical rather than (only) by physical parameters. A 
careful analys is is necessary to d i sentangle the respective influ- 
ence (see e.g..lPessah et al.ll2007HFromang & Papaloizoull2007t 



iRegev & Umurhanll2008l) r 

Global simulations, on the other hand, follow the evolution 
of the entire system, albeit with a much coarser resolution than 
local ones. Thus, they can account for the large-scale structure of 
stars and disks, for the back-reaction of the MRI instigated tur- 
bulence on the global flow, and allow one to draw conclusions on 
how the saturated state depends on global properties of the sys- 
tem, e.g., the density or pressure stratification. However, fore- 
going the ability to resolve short-wavelength modes, the growth 
of the MRI will be underestimated or suppressed even entirely. 
In many applications of numerical analysis, it is possible to use 
suitable models for the unresolved physics on sub-grid scales, 
e.g., sub-grid diffusivity. This requires a good knowledge of the 
physics on these scales, and is facilitated greatly if processes at 
the unresolved scales act merely as a sink for kinetic or mag- 
netic energy cascading down from the integral scale. During the 
growth of the MRI, however, the power shifts gradually from 
short to long modes. Thus, sub-grid models for global MRI sim- 
ulations tend to be complex, and are not used widely. 

As a remedy for this problem, global simulations may be 
performed using unrealistically strong initial fields to guaran- 
tee that the fastest growing MRI modes are resolved numeri- 
cally. This approach presumes that the unresolved MRI modes 
are able to amplify the much weaker actual initial fields to the 
field strengths used as initial value in global simulations. This 
assumption can be justified, if the MRI acting on the unre- 
solved scales saturates at the initial field strengths imposed in 



global simulations, i.e., if rapid amplification by the MRI takes 
place over many orders of magnitude. However, this can only 
be proven by high-resolution local simulations. Enhancing the 
initial magnetic field by a constant factor throughout the com- 
putational domain, as it is often done in global simulations, is 
problematic as the MRI is a local instability, i.e. it is not ex- 
pected to cause a constant amplification of the field everywhere. 
The ambiguities regarding differences between the topology of 
this artificially enhanced field and that of a field amplified lo- 
cally by the MRI add to the uncertainties clouding the influence 
of magnetic fields on the overall dynamics. 

Both the local and the global numerical approach has been 
used for studying the MRI in accretion disks, and this combined 
effort has led to the rapid development of the field. Simulations 
of the MRI in core collapse supernovae, on the other hand, have 
not yet reached this advanced stage, mainly because of the weak- 
ness of the initial field o f the progenitors. A ccording to current 
stellar evolution models dHeger et al.l 12005 ). the canonical pre- 
collapse magnetic fields are so weak that they are unable to affect 
the dynamics of the explosion unless they are amplified strongly. 
Correspondingly, the wavelengths of the fastest growing MRI 
modes are approximately a few meters at most Q Thus, the pos- 
sible importance of MHD effects in core collapse supernovae 
depends on the existence of mechanisms which can amplify the 
field efficiently during core collapse and the post-bounce phase. 
The timescale available for the growth of the magnetic field is 
set by the time required to turn the accretion of matter onto the 
proto-neutron star into an explosion, i.e., a few hundreds of mil- 
liseconds. As already mentioned above Akivamaetal. (2003|) 
suggested that the MRI might provide this mechanism. They 
estimated the saturation field strength to be 10'^ - 10'^ G, i.e., 
clearly in excess of the artificially enhanced initial field strengths 
used in global simulations. 

Up to now, there exist only global simulations of MHD core 
collapse supernovae, which evolve the entire core of a mas- 
sive star through gravitational collapse, bounce, and ex plosion 
(e.g., Kotake et al. 2004; Yamada & Sawai 2004; Takiwa kTetall 
[2004); lObergaulinger et al.l l2006b a; B urrows et al.ll2007l) . These 
global simulations fail to fin d the MRI unless the y employ dras- 
tically stronger initial fields. lObergaulingeret al . (2006b a), e.g., 
require a pre-collapse field strength exceeding 10'^ G to resolve 
the MRI in the post-bounce state. The rationale behind the arti- 
ficially increased initial field strengths is that, once triggered by 
the differential rotation in the proto-neutron star, the MRI will 
exponentially amplify a much weaker seed field up to the values 
used in the simulations. 

Due to the lack of local simulations, the importance of the 
MRI in MHD core collapse models remains unclear. As a first 
step to resolve this issue, we have performed high-resolution 
simulations of small parts of simplified post-bounce, rotating, 
magnetized cores. We have used a recently developed high- 
resolution MHD code, and employed s hearing-disk boundary 
conditions dKlahr & Bodenheimeiil2003h . These boundary con- 
ditions derive from the shearing-sheet boundary conditions of 
lHa wlev&Balbusldl992 ), but allow one to consider global gra- 
dients of, e.g., density or entropy. Combining elements of global 
and local simulations, viz. the presence of global background 
gradients, and a high-resolution local grid, we find it justified to 
call our approach semi-global (for more details see Sect. l3.2| ). 



' Note, however, that these predictions still involve uncertainties, and 
hence rare, but much more strongly magnetized progenitors cannot be 
excluded presently. 
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Differences in the physical conditions in disks and stars im- 
pede the direct application of the MRI results from accretion 
disks to supernovae. Most obviously, the geometry of both sys- 
tems differs strongly. Furthermore, while accretion discs are sta- 
bilized against gravity by (Keplerian) rotation, stars are sup- 
ported mainly by pressure gradients, with only a minor contribu- 
tion from rotation, i.e. thermal stratification is much more impor- 
tant in stars than in disks. Thus, entropy gradients can stabilize 
an MRI-unstable region or modify the instability in convectively 
unstable regions. Consequently, the problem of the MRI in core 
collapse supernovae has to be addressed by simulations account- 
ing for their specific properties, which is the goal of this study. 
We investigate the growth of the MRI from initial fields compa- 
rable to the ones expected from realistic stellar evolution mod- 
eling, and we seek to probe the possibility of MRI-driven field 
amplification under typical conditions of supernova cores and 
on timescales similar to the dynamic times of the system, (i.e., 
a few tens of milliseconds). Apart from the restrictions inher- 
ent to local and semi-global simulations, several simplifications 
limit our approach: we use simplified initial equilibrium mod- 
els, a simplified equation of state, and neglect neutrino heating 
and cooling. The main physical questions that we try to address 
are: (i) does the MRI grow on sufficiently short time scales to 
influence the explosion, i.e., within at most 100 msec, given typ- 
ical post-bounce rotation profiles and magnetic fields? (ii) How 
does an entropy gradient affect the growth of the instability? (iii) 
How does the saturated state of MRI-driven turbulence depend 
on these fa ctors? In particular, is the saturation field strength es- 
timated bv lAkivama et al.l (l2003b . i.e., the conversion of most of 
the rotational energy into magnetic energy, realistic? 

Analogous questions are studied by local simulations of the 
MRI in accretion disks. The answers may lead the way to for- 
mulate a turbulence model to be used in global simulations. The 
simplest model would provide a paramet rization of the angu- 
lar m omentum transport by an a viscosity ( Sha kura & Svunvaevi 
11973 ). i.e., a turbulent viscosity proportional to the local sound 
speed and the pressure scale height. However, despite a large 
number of local simulations, no unique formulation of an a 
model for accretion disks has been found up to now. Lacking 
similar comprehensive local simulations, a turbulence model for 
the MRI in supernovae is even less conceivable. Our simulations 
intent to provide only a first step towards these highly desired 
turbulence models. 

The paper is organized as follows: after a discussion of the 
main properties of the MRI in disks and stars (Sect.|2]i, we 
outline our numerical method in Sect.[3] discuss our results in 
Sect.m and summarize our main results and give conclusions in 
Sect.m 

2. MRI in discs and stars 

2.1. Physical model 

We work in the limit of ideal magnetohydrodynamics (MHD), 
solving the the equations of ideal MHD in the presence of an 
external gravitational potential tp, 

5,p + V,. [pv>] = 0, (1) 

d,p' + Vy [pV + Pj' - b'b'\ = pV,(/), (2) 

+ V,. [(e* + P*) - b'vtb'] = pv^V jip, (3) 
d,b =-cVxE, (4) 



Vjbj = 0. 



(5) 



Here, p, p, v, and e^. denote the mass density, momentum density, 
velocity, total-energy density, of the gas, respectively; b is the 
magnetic field. The total-energy density and the total pressure. 
Pi,, are composed of fluid and magnetic contributions: = e -n 
|pv^ + jh^ and P^, - P + jb^ with the internal-energy density 
E and the gas pressure P - Pip, s, . . .). The electric field, E, is 
given hy E - -^xb. Here, c = 2.998 x 10'" cm s"' is the speed 
of light in vacuum and Einstein's summation convention applies. 

We use the hybrid equation of state (EOS) due to iKeil et~al] 
(Il996h as a rough model for neutron-star matter. Following this 
EOS, the total gas pressure, P, consists of a barotropic part, P^, 
and a thermal part, Pth. The two parts are given by 



Pb 
Pth 



(Fth - l)£th- 



(6) 
(7) 



Here, Fb and k refer to the barotropic adiabatic index, and the 

polytropic constant ofthe EOS, respectively; fith = e-^'b/(rb-l) 
is the thermal part of the internal energy, and Fth the corre- 
sponding adiabatic index. Please note that we consider only sub- 
nuclear densities, p < pnuc = 2 x 10''* g cm"-^ here since the max- 
imum density reached in our models is a few times lO'"^ g cm"^. 
We define a pseudo-entropy S for this equation of state by 



(8) 



In the Schwarzschild criterion for convective stability (see be- 
low), this quantity appears instead of the entropy of, e.g., an ideal 
gas. 

A few quantities used frequently in the remainder of this pa- 
per are: 



1 . the Alfven velocity 
1*1 



CA = 



2. the (local) magnetic energy density 



Smag — 2 ' 



(9) 



(10) 



(11) 



and the corresponding volumetric mean value 

3. the (local) Maxwell Stress tensor 

Uij^bibj, (12) 
and the corresponding volumetric mean value 

Mu-\^ / ^^^^i- ^^^^ 

We will use most frequently the component M^t^ which 
governs the transport of angular momentum in radial direc- 
tion, and we will sometimes refer to this component as the 
Maxwell stress for short. 
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2.2. General properties of the MRI 

The stability c r iteria for the MRI w as first discovered by 
Velikhov' ( 1 959') : 'Chan drasekhad d 1 960 *) and further discussed by 
Balbus & Hawley ( in a series of papers. These authors an- 
alyze wave-like (WKB) perturbations of the form exp[i{k-r+a>t)] 
in a background equilibrium of the MHD equations. For con- 
venience cylindrical coordinates {m,4>,z) are used in the fol- 
lowing equations. From the dispersion relation, they derive the 
criteria for exponential growth and, if applicable, the growth 
rates of WKB modes. Because the main astrophysical context 
of this series of papers is accretion discs, some assumptions 
are made which considerably simplify the analysis: i) weak 
magnetic fields where |va| ^ min(cs, |v^|), ii) incompressible 
gas (Boussinesq approximation), and iii) angular velocity con- 
stant on cylinders, Q.(vj). The discussion is mostly restricted to 
thin discs (i.e., to equatorial regions and to a Keplerian rota- 
tion law) and 5/3-polytropes. Under these assumptions the sta- 
bihty criterion for a diff erentially rotating magnetized fluid is 
(lBalbus&Hawlevlll991h 



JL^ = md^n^ > 0. 



(14) 



If the criterion is not fulfilled only modes with (dimensionless) 
wavenumber 



k< 



(15) 



are unstable, where k = k ■ Va/^, Km = ^cnt ' ^a/^ and = 
Hj^IQ}. The (dimensional) grow th rate of the fastest growing 
mode is (iBalbus & Hawlevll 1 992h 



are the gravitational, rotational, and buoyancy terms, respec- 
tively, and Fi = SlnP/Slnpl ,. It is convenient for the mode 
analysis and if cos 9k +0Xo define the quantity 



C = 



cos^ 6k 

tan^ 6k - 2B^Q, tan 9k + Q^B^ + 'R^)in^ 



(24) 



where the curl of Eq. (|2]i, i.e. the vorticity equation, has been 
used to simplify the expression of C. Note that this quantity de- 
pends on the direction of the perturbation 9k, but not on k- itself. 
If cos 6k - 0, which corresponds to velocity perturbations par- 
allel to the rotation axis, the value of C diverges, but all length 
scales and growth rates of the discussion below are finite, and 
can be computed by taking the limit cos 9k — » 0. 

In the absence of a magnetic field, i.e. P - 0, the stabil- 
ity condition is simply C + 4 > 0, which is equivalent to the 
Solberg-H0iland stability criteria for a non-magnetized rotating 
fluid (Tassoul 1978). 

Because we want to study instabilities of magnetized fluids, 
we consider hereafter only the case k Q. Then the stability 
condition is C > 0, which corresponds to that of Balbus (1995), 



g S + K e^ > 
(gxe^)-(Sx'R) > 0. 



(25) 
(26) 



Modes with wave numbers smaller than the (dimensionless) crit- 
ical wavenumber 



<9k^PC 



(27) 



(16) 



which is independent of the magnetic field and corresponds to 
(dimensionless) wave numbers close to ^crit- 

However in the context of core collapse supernovae some of 
these assumptions do not apply: entropy and composition gradi- 
ents are important, more general rotation laws Q.{m,z) have to 
be considered, and the analysis can no longer be restricted to 
equatorial regions. In this general case t he dis persion relation of 
WKB modes is (cf. Balbus 1995; Urp iiill 19961) . 



{C?-k^f - [u 



■ P) + Gj\+A cos^ 



Ok) 



4p-cos^6k 







(17) 
(18) 



where d> = co/D. is the dimensionless growth rate of the insta- 
bility, and 9k is the angle between k and the z-axis. The (dimen- 
sionless) frequencies related to buoyancy terms and differential 
rotation are 



and 

-,2 



R 



1 
Q2 



^2 



k^ 



(k ■ e^)(k ■ K) 



P 



respectively, where e-^ is the unit vector in xu- direction, 

P 

Vp VP _ 1 51nP 
~p ~ ~ "Tx ds 



Vs 



(19) 
(20) 

(21) 
(22) 
(23) 



are unstable and grow. The critical wavenumber depends on the 
angle 6k in a complicated way involving, in general, the rotation 
profile, "R, the thermal structure, S, and the stratification, Q. 

Two branches of unstable modes arise from the dispersion 
relation with k Q (Urpin 1996.) : the branch of Alfven modes 
appearing for C < 0, and the branch of buoyant modes which 
only appear for C -H 4 < (Fig.[l]l. 

For a given 9k the fastest growing mode is obtained from the 
condition (J^w = 0. For -8 < C < it has a (dimensionless) 
wavenumber 



^FGM = cos 9k 



V=C(CT8) 



and a (dimensionless) growth rate 



(WFGM - COS ( 



(28) 



(29) 



If C < -8, the fastest growing mode corresponds to Pgm = 0, 
i.e. it is dominated by buoyant modes with a (dimensionless) 
growth rate 



WfGM - COS ( 



VC + 4. 



(30) 



Thus, there exist two different instability regimes depend- 
ing on the value of C. For -4 < C < only Alfven modes 
are possible. This mag neto-shear regime was discussed by 
iBalbus & Hawleyl (11991 ). A mixed regime is found for -8 <C < 
-4, where both Alfven and buoyant modes compete. For C < -8 
the buoyant modes completely dominate the growth of the insta- 
bility, and this regime is thus called magneto-convectivertgimt. 
It is similar to the convective regime, but the critical wavenum- 
ber is determined by the strength of the magnetic field. 
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beit with a vanishing growth rate as their wavelength A ap- 
proaches infinity. The growth rate peaks for 



-0.01 



Fig. 1. Imaginary part of the growth rate normalized to the imag- 
inary part of the maximum growth rate, 3(d>)/3(cjfCM) as a 
function of C and k I cos 0k = k ■ Va/C^^cosS^). The dashed 
line shows the value of k corresponding to the fastest growing 
mode, 3(cj)/3(cjfCM) - 1, the solid line gives the boundary be- 
tween the two branches of unstable modes (Alfven and Buoyant 
modes), and the dash-dotted line corresponds to the stability 
limit (k - kail). For -4 < C < only Alfven modes appear, with 
a narrow spectrum of fast growing modes close to k^ii (dash- 
dotted line). For C < -4 buoyant modes appear and become 
dominant for C < -8. In the latter case the spectrum of fast 
growing modes is much wider covering the entire region from 
^crit to 0. 



Note that for a given fluid element the behavior of the unsta- 
ble modes depends on the angle 0^. Thus, different regimes can 
hold in different directions. To find the absolute fastest growing 
mode of a fluid element, i.e. not considering a fixed angle 6^, one 
has to determine the zeros of dd>/d0k, which involves the solu- 
tion of a quartic equation. This fact makes a more detailed study 
of the instability difficult. 

To simplify the analysis, we restrict ourselves in the follow- 
ing discussion to regions near the equator, where it is reasonable 
to assume only a radial dependence of the hydrodynamic quan- 
tities and a vertical magnetic field. Therefore, 



(31) 



where = S-ffis the square of the Brunt- Vaisala or buoyancy 
frequency. Because C90 does not depend on 6^, all modes of the 
considered fluid element belong to the same branch of modes, 
i.e. they are either buoyant modes or Alfven modes. All modes 
with wavelengths shorter than 



'Icrit = 



2n 



2n |va| 



(32) 



are stabilized by magnetic tension. It is easy to show that the 
modes grow faster when k is parallel to the z-axis {0k - 0), 
i.e. velocity and magnetic field perturbations grow in direc- 
tion perpendicular to the rotation axis. The stability criterion, 
can easily be interpreted according to the relative 
size of the buoyancy term, A^^, and the shear term, 'R,;,. Several 
different regimes result (Fig.|2l): 

- Magneto-shear instability (MSI): <S A^^ and -4Q^ < 
A^^ + '??n7 < 0. All modes longer than /Icrit are unstable al- 



^MR\ = 27r/fcFGM ~ V2/lc: 



(33) 



where the limit |C| <K 8 is used to obtain the second ex- 
pression. It is important to note that /Imri, which scales with 
the background field strength bo, becomes small for weak 
initial fields. Hence, in the limit of a pure shear instability, 
only relatively strong initial fields are accessible by numer- 
ical simulations due to the restrictive constraint on the grid 
size imposed by the requirement to resolve /Imri by at least 
several grid zones. 

- Magneto-buoyant instability (MBI)-^N^ <«; iRj^ and -4Q^ < 

+'R„ < 0. This regime resembles the magneto-shear 
regime, but the instability is not driven by the shear, but 
rather by the unstable stratification. 

- Magneto-convective instability: <K "R^ and A^^ + "R^ < 
-4Q^. This regime corresponds to magnetized convective 
flow. The main difference is the stabilization of short modes 
(A < Adit) due to the magnetic tension. The more important 
the negative entropy gradient becomes with respect to the an- 
gular velocity gradient, the faster is the growth of infinitely 
long modes compared to the growth rate at /Imri- 

- Hydrodynamic shear instability: "Rcj ^ and A^^ -I- K^j < 
-4Q^. This case is not of interest in core collapse since for 
the differential rotation of PNS we always find > 1.5Q?. 

Core collapse occurs in general in a mixed regime, where Alfven 
modes and buoyant modes compete. Therefore, none of the 
above mentioned pure regimes holds for the MHD instabilities 
appearing during core collapse. 



3. Method 

3.1. Code 

We use a newly developed three-dimensional Eulerian MHD 
code (Obergaulinger et al., in preparation) to solve the MHD 
equations, Eqs. ([T]i-(|5]l. The code is based on a flux-conservative 
finite-volume formulation of the MHD equations and the 
constrained-trans port sc heme to maintain a divergence-free 
magnetic field (fEvans & HawleY|ll988h. Using high-resolution 
shock capturing methods (e.g., 'lLeVequel fT992h . it employs var- 
ious optional high-order reconstruction algorithms including a 
total-variation diminishing piecewise-linear (TVD-PL) recon- 
struction of second-order accuracy, a fifth-order monotonicity- 
preserving (MPS) scheme jSuresh & H uvnh 1997), and a fourth- 
order weighted e ssentially non-oscillatory (WEN04) scheme 
(iLevy et alJl2002h . and approximate Riem ann solvers ba sed on 
the multi-stage (MUSTA) method (Toro & TitarevI 12006). The 
simulations reported here are performed with the MPS scheme 
and a MUSTA solver based on the HLL Riemann solver dHartenl 
Il983h . 

The simulations are performed using cylindrical coordinates, 
and include both three-dimensional and two-dimensional (i.e., 
axisymmetric) models. The computational grid covers a region 
of a few (typically one or two) kilometers aside resolved by at 
least 26 and at most 800 zones per dimension, corresponding to 
a resolution between 40 and 0.62S m. 



- The reader should not con fuse this inst ability with the magnetic 
buoyancy or Parker instability iParkej[l966l) , related to the magnetic 
field strength gradients. 
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Fig. 2. Stability regions in the plane 11^/0^ vs N'/Q?. The 
sohd thick line separating the stable region from the magneto- 
rotational instabilities (MSI and MB I) corresponds to C = 0, and 
the solid thick line separating the magneto-rotational instabilities 
from the hydrodynamic instabilities (convection and shear insta- 
bility) corresponds to the C - -4. The mixed zone is arbitrarily 
defined by IK^/Q- - N^/Q?\ < 2. 
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Fig. 3. Hydrostatic structure of the initial models. The diagram 
shows the gravitational potential 1^19 - (p/{lQ^^ erg cm"^) (solid 
lines, right y axis), the density pn = p/(lO^^ g cm"-') (dashed 
line, left y axis), and the angular velocity Q.3 = f2/(10^ s"') 
(dash-dotted line, left y axis). The entropy profile of this specific 
model is assumed to be flat. 



3.2. Boundary conditions 

In local simulations, the choice of boundary conditions is a cru- 
cial issue, with possibly subtle effects on the flow geometry. The 
standard technique for local simulations of the MRI in ac cre- 
tion di sks is the shearing-sheet method due to lHawlev & BalbusI 
( II992I) . This approach consists of two important ingredients: (i) a 
transformation into a frame of reference co-rotating at the mean 
angular velocity of the shearing box, fJo, and the linearizion 
of the rotation profile around Qq; (ii) the use of shearing-sheet 
boundaries in the radial direction, and (in most cases) periodic 
boundary conditions in the perpendicular directions. 

Periodic boundary conditions are often used in simulations 
of a small, representative sub-volume of a larger system. These 
boundary conditions are based on the idea that the entire sys- 
tem is covered by a homogeneous (e.g., cubic) lattice of iden- 
tical sub-volumes. Consequently, the, e.g., left boundary of the 
simulated sub-volume is identified with the right boundary of an 
identical sub-volume translated by one lattice width. 

A shearing box represents only a small part of the entire 
system. The influence of larger scales is considered by suit- 
able boundary conditions, the most natural choice being periodic 
ones. These boundary conditions, however, do not allow one to 
impose global gradients throughout the shearing box, e.g., for 
differential rotation (dj^Q. + 0). This shortcoming is eliminated 
by the linearization of the rotation profile and the transformation 
into the co-rotating frame since, in this case, the deviation from 
the background profile, is the dynamical variable rather than 
Q itself. Thus, it is possible to use periodic boundary conditions 
in the radial direction accounting for differential rotation by an 
offset 6<p{t) - t {Q.oiit - Q/nn), as described by (Hawlev & Balbus 
Il992h . where are the angular velocities at the outer and 

inner radial surface of the shearing box, respectively. 



In contrast to accretion disks, thermodynamic variables in 
stars may have global gradients both in the direction per- 
pendicular and parallel to the gradient of O. Thus, standard 
shearing-sheet boundaries ca nnot be used. Instead, we follow 
iKlahr & Bodenheimed (l2003h and employ shearing-disc bound- 
ary conditions. We abandon the transformation into the co- 
rotating frame and assume radial periodicity of the deviation, 
6q = q- qo, of a variable q from a given background distribution 
qo, instead of periodicity of q itself. We define the background 
distribution qo by its distribution at the initial time f = 0, i.e. 
^0 = ?(w; t = 0). This recipe is applied to density, momen- 
tum, and entropy. As Klahr & Bodenheimer (2003), we observe 
the development of resonant radial oscillations which are sup- 
pressed, however, by damping the radial velocity in the first n 
(we use n - 2) computational zones at both radial boundaries. 
We point out that shearing-disc simulations allow for large-scale 
modifications of global gradients. In particular, angular momen- 
tum transport may modify the global rotation profile, and change 
the angular momentum and rotational energy of the matter in 
the computational volume. This process can eliminate the dif- 
ferential rotation causing the instability, and thus, terminate the 
growth of the MRI. 

As we will show later, the evolution of our models depends 
crucially on whether we do or do not apply this damping term. 
However, we note here that the artificial oscillations prevented 
by the damping do not have a strong influence on the evolu- 
tion of the MRI. We use, if at all, a damping of v-^ by 1.25% 
in the innermost and outermost two zones of the grid, which 
is a considerably weaker damping than in the simulations of 
Klahr & Bodenheimer (2003). Despite its relative weakness, the 
damping term is able to suppress weak radial motions across 
the boundary. Thus, it introduces a preferred length scale (the 
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radial size of the box) into the otherwise shear-periodic simula- 
tion. Comparing simulations with and without damping (we will 
refer to these boundary conditions by d and p, respectively), we 
can study the influence of a preferred scale on the MRI. 

The box size of standard shearing-sheet simulations does not 
define a preferred length scale, i.e., these simulations are scale 
free and entirely local. In shearing-disc simulations, in contrast, 
the scale height of the thermodynamic variables introduces a 
physical length scale into the simulation. If this preferred length 
scale is smaller than the entire size of the star or disk, the simu- 
lations can be characterized as being semi-global. 

The semi-global approach falls in between a purely local and 
a global one sharing merits and drawbacks with both methods. 
Similar to local simulations, semi-global simulations allow one 
to resolve a small part of the entire system better. Because they 
rely on a fixed lab frame and do not eliminate the mean rotation, 
the basic time scales are the same as in a global simulation of 
the same resolution. In a Keplerian disk dominated by rotation 
this might add a major difficulty to the numerical treatment of 
the problem. On the other hand, with pressure dominating over 
kinetic energy, the time step of our simulations is governed by 
the sound speed rather than the rotational velocity. As there is no 
way of eliminating the sound speed, we do not feel a need to use 
a shearing-sheet transformation. 

We expect the MRI in core collapse to grow and reach 
saturation within several tens of milliseconds. The time step 
6t < 6x/cs, on the other hand, is much smaller because of the 
high value of the sound speed in a post-collapse core (cs ~ 
lO'^cms '), where 6x is the width of the computational zone. 
Thus, we have to perform a large number (typically several mil- 
lions) of time steps, which implies a limit on the grid resolution 
we can afford in the simulations, although the resolution is still 
much better than that of a global simulation. 

3.3. Initial conditions 

We use equilibrium initia l model s based on post-bounce cores 
fromlObergauling er et al.l (I2006H) . Several tens of milliseconds 
after the core bounce, the shock wave has reached distances of 
a few hundred kilometers, the post-shock region exhibiting a 
series of damped oscillations as the proto-neutron star relaxes 
into a nearly hydrostatic configuration. We extract the radial 
profile of the gravitational potential along the equator of their 
model A1B3G3Q, and construct from that the density stratifica- 
tion within our shearing box solving the equation of hydrostatic 
equilibrium 



- pd^(p - d^P + vjpQ. 
for a given rotation profile 
□(tir) = QoZEr"" , 



(34) 



(35) 



where Q.o and oq are constants. The pressure is determined using 
the hybrid equation of state in the form 



assuming an entropy profile of the form 
S (m) - S + S livj - vjq) 



(36) 



(37) 



with constants S o and 5 1 . Eq. ( l34l ) is solved in a radial domain of 
size, AtiT, which is either one or two kilometers large, centered 



at vjo = 15.5 km. The structure of an initial model, characterized 
by the set of parameters flo - 1900 s"', = -1 .25, 5o = 0, and 
5 1 = 0, is shown in Fig. [3] This model has a radial density scale 
height of //p - -rf-p ~ 3.8 km, i.e., our computational grid covers 
a significant fraction of a density scale height. The rotation rate 
of ~ 2 000 s"' corresponds to that of a rapidly rotating proto- 
neutron star with a rotational period of ~ 3 ms. 

Assuming that the background gravitational potential is a 
function of m only, we construct cylindrically symmetric ini- 
tial models. This approximation is justified by the small size of 
the simulation box in z-direction (1 km) compared to its radial 
position (15 km). 

We added three different types of initial magnetic fields to 
the initial hydrostatic model: 

Model UZ: a uniform B-field in z-direction, b — (0, 0, /Jq) . 
Model VZ: a B-field in z direction with vanishing net flux, b = 
(0, 0, sin (27r(m - mo)/Am)f. 

In all models, the initial field is weak, both in compaison with the 
thermal and the rotational energy of the models. The weakness 
of the associated Lorentz force justifies the use of hydrostatic 
instead of hydromagnetic equilibria as initial conditions. 

From Eq. ( |33] ) and the values of typical model parameters we 
obtain the following estimate for the MRI wavelength 



/ipGM ~ 6.9 km 



P 



10'5G/\2.5 lO'^^gcm- 



1900 s- 



■(38) 



To properly simulate the evolution of the MRI, /Imri should be 
resolved by at least a few grid zones. Using a grid resolution 
of 10 m or 20 m, we thus can follow the gr owth of the MRI for 
magne tic fields exceeding several 10'^ G dObergauUnger et af] 
I2006bllal) . To trigger the instability, we impose a small random 
radial velocity perturbation with an amplitude of a few times 
2Q-4..-2 q£ j.jjg rotational velocity. 



4. Results 

4.1. General considerations 

In axisymmetry, the growth of the MRI requires a non-vanishing 
poloidal initial field. Axisymmetry restricts the dynamics of the 
MRI, suppressing a class of instabilities that affect the evolu- 
tion of MRI-unstable modes (see below). Consequently, the pre- 
dictive power of axisymmetric simulations for the evolution of 
the MRI is limited, and we cannot rely on them in determining 
the saturation amplitude of the instability in supernova cores. 
The growth of the instability does, however, not differ strongly 
from full 3D models. Thus, we can use 2D models to determine 
growth rates, while detailed conclusions can only be drawn from 
3D models. 

In axisymmetry, the flow is dominated by channel modes, 
a pattern of predominantly radial flows of alternating direc- 
tion stacked in z direction (Balbus & Hawlev 1991). As the MRI 
grows, the channels start to merge and their characteristic length 
scales increase, but they survive as coherent flow structures 
throughout the entire evolution and, particularly, do not dissolve 
into turbulence. 

The analysis of iGoodman & Xul d 19941) shows that channel 
modes are an exact nonlinear solution of the axisymmetric MHD 



^ This model experiences a core collapse that is halted by the stiffen- 
ing of the equation of state above nuclear matter density. 



In general, the channels are oriented parallel to the gradient of Q., 
wherever it points to. 
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equations, which explains their stabihty observed in many nu- 
merical simulations. They are, on the other hand, unstable to 
genuinely 3D parasitic instabilities of, e.g., Kelvin-Helmholtz 
type. Consequently, in 3D, the channel modes appearing during 
the early growth phase of the MRI, do not persist until satura- 
tion. Instead, the channels decay due to the growing parasitic 
instabilities, and turbulence develops. 

This basic picture emerged from many simulations of the 
MRI in accretion disks. As we will discuss in the following, our 
simulations confirm this result for the MRI in supernova cores. 

4.2. Axisymmetric models with no entropy gradient 
4.2.1. Uniform initial magnetic fields 

Our models having no entropy gradient show the same dynam- 
ics as that observed in previous simulations of the MRI in accre- 
tion discs (see, e.g., Balbus & Hawlev 1998). We discuss first the 
models with a uniform initial field in z-direction (model series 
UZ2) focusing on models with a rotational law given by Qq - 
1900 s"' and an = -1.25 (see Eq|35]l. The evolution of these 
models is characterized by an exponential growth of the mag- 
netic field, see e.g., Fig.|4]for a model with bf^-2x lO'^ G. The 
fastest growing MRI mode is well resolved in this model, and 
its growth rate ctmri = 1-08 ms"' is found to be close to the the- 
oretical prediction crjyiRi = 3{a>fCM) ~ |Q'nfio/2| - 1.14ms"' 
(see Eq. (|29]l). The magnetic field reaches a maximum value of 
about lO'^G at f a; 15 ms, and the mean Maxwell stress com- 
ponent Mo70 (see Eq. (fT3] )) becomes large enough to alter the 
rotation profile considerably within a few tens of milliseconds. 
Consequently, the angular momentum of the gas drops drasti- 
cally at f a 25 ms. 
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Fig. 4. Evolution of the mean magnetic energy density e™^ 
(solid black line), the mean energy densities corresponding to 
the m (dotted red), (dashed brown), and z (dash-dotted green) 
component of the magnetic field, and the absolute value of the 
mean Maxwell stress component M^,p (dashed blue line) for 
an axisymmetric model with an initially uniform magnetic field 
= 2xlO'^Gin z-direction, and a rotation law given by 
Qo = 1900 s"' and an - -1.25. The model was computed in 
a box of L^x - 1km x 1km with a grid resolution of 5m. 



The growth of the MRI proceeds via channel flows, whose 
vertical extent and number depends on the initial magnetic field. 
Two typical channel flows are shown in Fig.|5] During the early 
phase of the instability {t - 10.6 ms; upper panel) eight distinct 
channels are present each one consisting of a pair of up- and 
down-flows in radial direction. The magnetic field is organized 
into eight elongated radial sheets, and this pattern is also im- 
printed onto the distribution of Q, as the magnetic field enforces 
co-rotation along field lines. 

A flow topology dominated by channel modes implies a 
phase of exponential growth of the magnetic field, which ends 
when the channel modes are disrupted and a less organized, more 
turbulent state ensues (in Fig. |4] this happens at f » 11 ms). In 
most axisymmetric models, the turbulent state is only of tran- 
sient nature, because after some time coherent channel flows 
form again leading to a secondary phase of exponential growth 
(see Fig.|4]at t ^ 23 ms). 

Very late in the evolution (f = 30.7 ms; lower panel of 
Fig. Is) we find only one large-scale channel flow which ex- 
tends across the entire domain in radial direction. The mag- 
netic field is now predominantly radial and is concentrated near 
the channel boundary. This coherent flow pattern is the result 
of a strong transport of mean angular momentum by Maxwell 
stresses. The stresses enforce co-rotation along field lines, and 
consequently turn the rotation profile, initially constant on cylin- 
ders m - const., by 90 degrees, so that Q becomes a function of 
z only. We can distinguish two regions of slow and fast rota- 
tion inside and outside z e [-0.15; 0.25], respectively. Inside the 
slowly rotating channel matter is accreting towards the center 
with Vjjf ~ 4 X 10^ cm"', while the rapidly rotating gas outside 
the channel has much slower, random velocities. 

To investigate the dependence of the channel geometry on 
the initial magnetic field and the grid resolution we compute 
Fourier spectra of the radial component of the magnetic field, 
bj^, for models with Qq - 1900s"' and an = -1.25, and an 
initial field strength of 4, 10, and 20 x 10'^ G, respectively. The 
simulations are performed in a box of either 1x1 km^ or 2x2 km^ 
using a 400^ grid (Fig.|6]l. At each radius we Fourier-transform 
btj(z), and the resulting spectra bj^{k~) (where is the vertical 
wave number) are then averaged over radius. We applied this 
procedure to the models during the growth phase of the instabil- 
ity at f a! 7.5 ms. The first set of models (1 km^ domain, 2.5 m 
spatial resolution) exhibits growth rates close to the theoretical 
values, while this only partially holds for the models of the sec- 
ond set of models (4km^ domain, 5 m spatial resolution). Due 
to insufficient spatial resolution the MRI in the model with the 
weakest initial field {b^ = 4 x 10'^ G) grows slower than theoret- 
ically predicted. However, for the two more strongly magnetized 
models of this set {b^^ - 10, and 2 x lO'-' G) the fastest growing 
modes are well resolved, and the MRI growth rates agree with 
the theoretical ones. 

For each model the spectrum shows a distinct maximum cor- 
responding to a dominant vertical length scale given by the width 
of one channel mode. The position of this maximum is a func- 
tion of the initial magnetic field only, A;inax ^o ' ^^^^ dots, 
neither depend on the size of the computational domain nor on 
the resolution. A dependence on the last quantities is only ob- 
served, if the fastest growing mode is under-resolved. In this 
case, we recover the low-A: wing of the spectral peak, but find 
a truncated spectral distribution at higher wave numbers/smaller 
length scales. 

MRI theory predicts that the growth rate is independent of 
the initial field strength. Neglecting magneto-convective modes, 
we can expect to observe this behavior in numerical simulations 
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Fig. 5. The channel modes present in two snapshots taken from 
the model for which Fig. |4] shows the time evolution. The snap- 
shots are taken at f = 10.6 ms (upper panel) and t - 30.7 ms 
(lower panel), respectively. The panels show the color coded an- 
gular velocity Q, the magnetic field lines (white), and the flow 
field. The colors of the velocity vectors indicate the magnitude 
and the direction of the flow: up- and down-flows are represented 
by blue and red vectors, respectively, their color intensity corre- 
sponding to the absolute value of the (poloiodal) velocity (the 
darker the larger). The maximum velocities are 2.7 x 10^ cm s"^ 
(upper panel), and 8.8 x 10** cms"' (lower panel), respectively. 
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Fig. 6. Radially averaged Fourier spectra of the radial component 
of the magnetic field, b-^iz), for diff'erent uniform-field models 
at t X 7.5 ms. Models with /jq = 4 x lO'^G, lO'^G, and 2 x 
10'^ G are shown by the red dash-dotted, the blue dashed, and 
the black solid line, respectively. Thick and thin lines refer to a 
computational domain of 1 km^ and 2 km^, respectively. For all 
models a grid of 400^ zones is used. 



growth rate is smaller than crpcM, because the MRI wavelength, 
i.e., the wavelength of the fastest growing mode, exceeds the box 
size. Thus, we can only properly simulate the slower growth of 
shorter modes. 
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only if the grid is sufficiently fine to resolve the fastest grow- 
ing modes close to /Imri- Otherwise, if the grid is too coarse the 
growth rate should be much smaller Our simulations reproduce 
this behavior We show a comparison of the maximum growth 
rates from Unear analysis (o"fgm = 3('j^fgm) ~ Ijffw^^ol) and 
the numerical ones for models with different initial rotation laws 
(f2o ranging from 950s"' to 1900s"', and an from -1 to -1.25) 
in Fig.|7] If /Imri is under resolved for a given initial field ^o, the 
growth rate increases with ba, but once the MRI wavelength is 
well resolved, the growth rate becomes constant as theoretically 
predicted. Fig. |7] implies the following criterion for a sufficient 
resolution of the MRI: Atu > 2D.o/ca- The growth rate of the 
instability does not depend on the size of the computational do- 
main. For models with strong initial fields the computed MRI 



0.0 I i , 

1 10 
Ca / ( 8 X ) 

Fig. 7. Growth rate cr of the MRI for axisymmetric mod- 
els with uniform initial field as a function of the initial 
Alfven speed normalized to the rotational velocity and the 
grid resolution, ca/(5xQo). The colored symbols distinguish 
different initial rotational laws, where (Q.{),an) are equal to 
(1900s-', -1.25) [black plus signs], (1900 s"',-!) [red di- 
amonds], (950 s"', -1.25) [green squares], and (950s"',-l) 
[brown triangles], respectively. 



10 



M. Obergaulinger et al.: MRI in core collapse supernovae 



4.2.2. Channel disruption and IVIRI termination 

As long as the dynamics of the model is dominated by channel 
modes, the MRI grows exponentially. We observe a termination 
of its initial exponential growth - henceforth called MRI termi- 
nation - as soon as the coherent channels are disrupted. Further 
MRI growth occurs after an eventual reformation of the channel 
flows. To understand these processes better, we study MRI termi- 
nation in a large number of axisymmetric models with different 
initial magnetic fields, boxes of different size and grid resolution, 
and diff'erent boundary conditions. 

Fig. [8] shows the value of the mean Maxwell stress compo- 
nent M"^™ at MRI termination as a function of the initial mag- 
netic field strength, bo for models with a rotational law Q - 
1900 s^'tu"'^^ and a vanishing entropy gradient. We can dis- 
tinguish two classes of models according to the boundary condi- 
tions applied in the simulations (see section [T2] l, the qualitative 
difference between the models with and without velocity damp- 
ing near the boundaries being quite remarkable given the weak 
damping we apply. In models with damping M^™ grows with 
increasing initial field strength until it levels off at a grid size 
dependent value (colored bands in Fig. [8]). For the same models, 
when simulated without damping, we find that M^™ oc b]^^^'' 
(gray band in Fig. [8]) independent of the grid size. The "outliers" 
in the upper part of the figure correspond to models computed 
with a higher grid resolution than most other models. We will 
discuss this fact below. 

To determine whether the radial or the vertical size of the 
computational grid is responsible for the leveling off of the 
Maxwell stress in the runs with radial damping we simulate two 
models with a grid of 0.5 kmx2 km (for short called high models 
in the following), and 2 km x 0.5 km {long models), respectively. 
Our results show that the determining factor for the growth is pri- 
marily the radial rather than the vertical box size, as both models 
follow the behavior of M^.^ as a function of b^ for the respec- 
tive radial grid sizes. The two classes of models also exhibit re- 
markably different post-growth dynamics. In the high models, a 
few channel modes reappear from the turbulent state, and a sec- 
ondary phase of exponential growth of M^,,^ sets in. Eventually, 
two of the newly formed channel modes merge. By this process, 
which occurs repeatedly, the number of channels decreases, and 
the final state of the flow is dominated by one short but wide 
channel mode. In the long models, on the other hand, no sec- 
ondary exponential growth is observed, and the Maxwell stress 
remains approximately constant, albeit oscillating considerably 
due to the temporary presence of coherent flow patterns. 

To interpret these results, one has to analyze the mechanism 
responsible for the disruption of the channel modes. We discuss 
this mechanism for an undamped model with f2o = 1900s ', 
an = - 1.25, and an initial magnetic field strength feo =4xlO'-'G 
using a box of 0.5 km x 0.5 km and a resolution of 100 x 100 
zones. During the growth of the instability a few large chan- 
nels are present, which are disrupted at MRI termination (at 
* 15.9 ms). 

Fig. m illustrates the disruption of one of the channel flows 
in some detail. At f = 15.850 ms, the channel flow is still in- 
tact (left panel), and one recognizes two broad streams of in- 
flowing and out-flowing gas both permeated by a strong radial 
magnetic field of opposite polarity. A broad current sheet sep- 
arates the two flow regions. Owing to small-scale fluctuations 
in the flow, the field lines are not perfectly (anti-)parallel, and 
the current sheet is slightly deformed. These deformations act 
as seed perturbations for resistive instabilities of the tearing- 
mode type. Although we evolve the equations of ideal MHD 
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Fig. 8. Volume-averaged Maxwell stress component at 
MRI termination as a function of the initial magnetic field 
strength, bo, for axisymmetric models with uniform initial 
magnetic field in z-direction , a rotational profile O = 
1900 s ' w '-^^, and a vanishing entropy gradient for a set of 
axisymmetric models. Blue, green, orange, and red symbols cor- 
respond to models computed in a square box having an edge size 
of 0.5, 1, 2, and 4 km, respectively. Models computed with and 
without velocity damping at the radial boundaries are denoted by 
asterisks and diamonds, respectively. The latter models show a 
box size independent scaling M^™ ~ (b^y^^'' (gray band), while 
in models with damping M'^™ saturates at high field strengths 
the saturation value depending on the box size (colored horizon- 
tal bands). 



neglecting resistivity, the presence of numerical resistivity en- 
ables the growth of these instabilities, leading to a reconnection 
of anti-parallel field lines. As a consequence, the elongated cur- 
rent sheet dissolves into a configuration of X and O points (lo- 
cated at tir ~ 15.25 km and ~ 15.5 km, respectively; see middle 
panel of Fig.|9]l. When field lines reconnect near the X point, the 
fluid is accelerated away from the reconnection point towards the 
O point. This causes the intense gas flow in positive radial direc- 
tion at (ziT,z) ~ (15.35, 0.1)km. The change of the topology of 
the magnetic field and of the flow continues shortly afterward 
(at t - 16.078 ms; right panel of Fig.|9]l. The O point has grown 
in size, and the fluid is in vortical motion. As the vortex grows, 
field lines in the vortex are advected towards field lines of oppo- 
site polarity belonging to an adjacent channel flow (centered at 
z ~ -0.15 km), where reconnection occurs. Note the formation 
of a second X point at (ct,z) - (15.38, -0.03)km (Fig.|9l right 
panel). 

To demonstrate the growth of the tearing-mode instability 
and to support its importance for MRI termination, we compare 
the evolution of the mean Maxwell stress component M^jj^ and of 
the magnetic energy density of the z-component of the magnetic 
field, e^^„ of the model (see Fig.fTOli. Before the tearing mode 
grows M^^ and e^^^ increase, their growth rates being similar 
to that of the MRI. At t - 15.850 ms, the growth rate of ej^^^g 
becomes larger than that of the MRI by one order of magni- 
tude within less than 0.2 ms, whereas M^j^ approaches a maxi- 
mum. Once the tearing mode is fully operative (f = 15.983 ms), 
the growth of ejjj^g becomes slower but still continues due to the 
appearance of more tearing modes (see, e.g., the right panel of 
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Fig. 9. The disruption of a channel mode in an axisymmetric uniform-field model. We show a section of a model with an initial field 
bo = 4x 10'^ G computed on a grid of 0.5 x 0.5 km^. The left, middle, and right panels display the color-coded radial component 
of the magnetic field (top) and the current density = (V x b)^ (bottom) before {t = 15.850 ms), during (f - 15.983 ms), and 
after (t - 16.078 ms) the violent disruption of the channel flow, respectively. Additionally, magnetic field lines (black lines), and the 
velocity field (arrows; top only) are shown. The arrows are color-coded according to the magnitude and direction of the flow. Inflows 
and outflows are shown by gray and green vectors, respectively. The longest vector corresponds to a velocity of |v| = 7 x 10*^ cm/s. 



Fig.[TO]at (TiT,z) ^ (15.38, -0.03) km). Subsequently, e^3g begins 
to decrease as the tearing modes saturate. 
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Fig. 10. Temporal evolution of the absolute value of the mean 
Maxwell stress component Mj^^ (solid line; the line is colored 
black where M^^ < 0, and green otherwise), and of the mag- 
netic energy density of the z-component of the magnetic field, 
e;^2„ (dashed red line) of the model shown in Fig.|9l The vertical 
yellow lines mark the times of the snapshots shown in Fig.|9] 



According to the previous discussion the dynamics of the 
channel flows is dominated by the interplay between their 
growth due to the MRI and their destruction by resistive insta- 
bilities. Channel flows are unstable against tearing-mode-type 
instabilities at any point in their evolution. To study these insta- 



bilities in more detail, we have performed a set of simulations 
(see App.lAl) using simplified models of channel flows. We reca- 
pitulate our results, summarized in Eq. ( lA.llt , here: 



(39) 



where cTj, ca, cs, a, and 6x are the growth rate of the instabil- 
ity, the Alfven velocity corresponding to the channel magnetic 
field, the sound speed, the width of the channel, and the grid 
resolution, respectively, because the instability is not based on a 
physical resistivity, but is of purely numerical origin, no physi- 
cal transport coefficient appears in Eq. ( [39] l. However, our results 
can be interpreted in terms of an effective resistivity c^Sx, as de- 
tailed in App.lAl In our models the width of the channel flow, a, 
is set by the MRI wavelength corresponding to the initial vertical 
magnetic field: 

~ (40) 



^MRi ^ by ^J'R^ oc by{ ^/a^n^)) 

(see Eq. (l32l l. Eq. (l33T l, and Eq. (l35ll). The width remains con- 
stant during the growth, as only mergers of adjacent channels 
occurring as a result of resistive instabilities can change the field 
topology. 

Our basic proposition for MRI termination is that channel 
flows are disrupted once the growth rate of the resistive instabil- 
ity exceeds the MRI growth rate: 



o"r > cmri => MRI termination . 



(41) 



Using in addition the functional dependence of 0", (Eq. (39[). 
we can establish scaling laws for MRI termination for a given 
hydrodynamic background model. As the channel width scales 
with the MRI wavelength, a oc Amri, and as the MRI growth rate 
is given by itmri = 3(cjfCM) an^o (see Eq. (|29]l), we find for 
the Alfven speed at MRI termination 



(csf''{bf\nor''H6x) 



-4/7 



(42) 
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Fig. 11. Average Alfven velocity - of models with uniform initial magnetic field and without velocity damping near the radial 
boundaries - corresponding to the radial magnetic field at MRI termination normalized to {dxy'^^'' (left) and (h})^^'' (right) as a 
function of the initial magnetic field strength, (left), and the grid resolution 5x (right). The dashed lines represent the power 
laws expected from our analysis of resistive instabilities (see Eq. (l42b ). Note that we only consider well-resolved models (ctmri > 
0.95 ms-i) here. 



and for the corresponding Maxwell stress 



Ml 



(43) 



The latter equation implies that M^™ decreases with faster rota- 
tion. Two effects play a role in explaining this behavior. Firstly, 
slower rotation leads to slower MRI growth (ctmri ^ fio), and 
hence weaker magnetic fields are required for the tearing modes 
to overcome the MRI growth. Secondly, slower rotation implies 
wider channel flows (a oc (Qq)"'), i.e., resistive instabilities grow 
slower as cr,. oc a^^ oc (Qq)^ (see Eq. ([39]l). 

The qualitative features of these scaling relations are: 

1. Stronger initial vertical fields and correspondingly 
wider channels tend to suppress resistive instabilities. 
Consequently, MRI termination requires more strongly 
magnetized channel flows. 

2. Finer grid resolution implies less numerical viscosity, and 
hence larger values for cf ™ and M"'™. 

3. The scaling of M^™ with the sound speed implies a pro- 
portionality between the Maxwell stress and the background 
pressure: P oc Cg, and thus M|^™ oc f ^Z^. This scaling is rem- 
iniscent of the a-law in accretion discs according to which 
the (MRI-generated) viscosity is proportional to the gas pres- 
sure. 

Fig. [TT] shows that the average Alfven velocity correspond- 
ing to the radialfl magnetic field at MRI termination is well de- 
scribed by the scaling law given in Eq. ( l42b . Similarly, Fig.[8]and 
Fig. [12] (upper panel) confirm that the data for M^™ obey the 
corresponding scaling law (Eq. (|43])), too. The upper panel of 
the latter figure shows M^™ as a function of the initial magnetic 
field strength for models with different initial rotational laws. 

/ .\16/7 

Obviously, the proportionality M^™ k ib^j (light gray lines) 



' The radial field is typical for all three components. Thus the Alfven 
velocity corresponding to the total magnetic field shows the same de- 
pendence. 



provides a good approximation to the behavior of the models. 
Due to the small number of models the results should be taken 
with care, but a strong anti-correlation of Mj^™ with Oq is sug- 
gested. The data also do not support any dependence of M^™ 
on an- Finally, in Fig.[8]we noticed earlier some outliers at large 
values of Mj^™ which correspond to models computed on a fine 
grid. However, considering that MRI termination depends on 
grid resolution, all models lie within a narrow band which cor- 
roborates our scaling laws and provides more evidence of the 
importance of resistive instabilities in understanding the MRI. 
Consequently, physical (instead of numerical) transport coefli- 
cients should be used in MRI simulations, which may give rise 
to different scaling laws considering the growth rate of tearing 
modes. 



Models with velocity damping In models with uniform initial 
magnetic fields where the radial velocity is damped near the in- 
ner and outer radial boundary, i.e., in models located in the hor- 
izontal bands in Fig. [8] MRI termination happens earlier than 
predicted by our scaling laws, and the Maxwell stresses saturate 
for strong initial magnetic fields, the saturation value of M^™ 
being smaller for slower rotation (see lower panel of Fig.[T2li. 
This is due to the reconnection instability occurring close to the 
radial boundaries well before the theoretically predicted time 
of MRI termination. This premature reconnection is caused by 
the field geometry: due to the suppressed motion across the in- 
ner and outer radial boundary field lines must bent there in z- 
direction. Consequently, field lines of opposite polarity approach 
each other much earlier than in models without velocity damp- 
ing, and eflicient reconnection ensues. In this case, the onset of 
reconnection is determined by the field geometry rather than by 
the initial field strength. With reconnection occurring in the bent 
flux sheets near the radial boundaries instead between parallel 
flow sheets in the bulk volume as for non-damping boundaries, 
the width of a flux sheet is less important in determining the 
resistive growth rates. Thus, a slower MRI growth and smaller 
Maxwell stresses are found for slower rotating models when ve- 
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Fig. 12. IVIaxwell stress, M"^™ at [VIRI termination as a function 
of the initial magnetic field strength, If^ for models with differ- 
ent initial rotational profiles. The upper and lower panels show 
models with non-damping and damping boundary condition, re- 
spectively. The colored symbols distinguish different initial ro- 
tational laws where (Qq, an) are equal to (1900 s -1.25) [blue 
diamonds], (1900 s"',-!) [green diamonds], (950 s"', -1.25) 
[brown asterisks], and (950 s"', -1) [red asterisks], respectively. 
Note that only models with a box size x = 1 x 1 km^ and a 
resolution of 50, 100, or 200 zones (per dimension) are consid- 
ered here. The light gray lines in the upper panel illustrate power 



laws oc [h~^ 



locity damping is imposed. Apart from the dependence of M^™ 
on Qo, we also find a dependence on an. Both dependences to- 
gether give rise to a monotone relation between the strong-field 
limit of M^™ and the MRI growth rate, ctmri, which qualita- 
tively agrees with the above reasoning. 

Once the initial channel flows are disrupted and the field ge- 
ometry is changed by reconnection, the mean magnetic and ki- 
netic energies, and the absolute value of the Maxwell stresses 
begin to fluctuate strongly around roughly constant values (see 
the phase between 1 1 ms and 23 ms in Fig.|4]i. Subsequently, a 
second phase of exponential MRI growth is possible, exhibiting 
a similar dynamics but involving less channel flows than the pre- 
vious growth phase. The reduced number of channels is probably 
due to the strong increase in the vertical magnetic field during 
the growth of the tearing modes. Similarly to their predecessors. 



the newly formed channels are also unstable against resistive in- 
stabilities, but due to their larger width their disruption requires 
much higher Alfven velocities, i.e., the MRI can lead to much 
higher Maxwell stresses in the second generation of channels. 
In principle, this process of formation and merger of channels 
can continue until only one single channel flow remains cover- 
ing the entire box. We note that in later growth phases, the radial 
velocity and the magnetic field strength are typically so large 
that damping at the radial boundaries, if applied, does not lead 
to early saturation. 

£1 [ s-M 



1736 1792 1848 1904 
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Fig. 13. An early state (f = 12.1 ms) in the evolution of a model 
with i7Q = 2xlO'-'G computed in a box of L„ x L, = 1 km^ 
covered by 200 x 200 grid zones. Shown are the same variables 
as in Fig.|5] The maximum velocity is 1 . 1 x lO^cm s ' . 



4.2.3. Models with non-uniform initial magnetic fields 

Models having a non-unifo rm initial magnetic field exhibit a 
different evolution (see also iBalbus & Hawlevlll998l) . To study 
this evolution we simulated a set of models varying the initial 
magnetic field configuration and the boundary condition (apply- 
ing velocity damping or not; see previous subsection). All mod- 
els rotate initially according to the law given in Eq. dSST l with 
Qo = 1900 s-' and an = -1-25. 

We considered three types of non-uniform initial magnetic 
fields all of which have only a z-component. The first one 



^lNF = ^()Sin(^)x 



(44) 



varies sinusoidally with radius and scales in addition as oc to 
guarantee that the net magnetic flux through the surfaces of the 
computational box at z = Zo,i vanishes. This is the standard zero 
flux field used in most MRI simulations. The second type of a 
non-uniform initial field considered by us is given by 



'^Abs rZNFl 



(45) 



Finally, the third type also has a vanishing net magnetic flux, as 
bi^„, but a step-like dependence on m, i.e.. 



^step = ^0® - ^c) — : 



(46) 
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where and tUc denote the Heaviside step function and the ra- 
dial coordinate of the center of the box, respectively. 

For the first type of models the MRI starts growing via a 
multitude of channel modes giving rise to less amplification of 
the magnetic field than in models with a uniform initial field. 
Separate channels develop in the two radial regions of negative 
and positive (see Fig.fTSl). whereas the channels span the full 
radial extent if the initial field is uniform. The channel flows 
do not merge to form a few large-scale channels, but are de- 
stroyed by turbulence. After reaching a transient maximum, the 
magnetic energy and the Maxwell stress level off at values much 
less than for uniform-field models. The magnetic field becomes 
strongest right after MRI termination (~ 10'^ G). After 60 ms the 
maximum field strengths are about 2 x lO'^^G, and decrease to 
10'"^ G until the end of the simulation. Fields of this strength can 
change the O profile significantly only on time scales of many 
tens of milliseconds, i.e., at the end of the simulation the rotation 
profile is basically unchanged. 

Fig.[T4lshows a comparison of the maximum Maxwell stress 
at MRI termination for the models with a non-uniform initial 
magnetic field. When velocity damping is applied the value of 
M^™ does not depend on the initial field geometry, and is deter- 
mined by reconnection of anti -parallel field lines occurring close 
to the boundaries. If no velocity damping is applied, the evolu- 
tion is similar to that of uniform field models, but depends on the 
field geometry. This finding can be understood in the light of our 
previous discussion of re-connective instabilities, and by the fact 
that M"^™ is determined by reconnection in the bulk volume, and 
not by reconnection near the boundaries. 

Models without velocity damping and the initial field, b~j^^^ 
develop large Maxwell stresses which increase with the initial 
field strength. The evolution of these models and the geometry of 
their channel flows are similar to those of models with uniform 
initial fields, i.e., the growth rates of tearing-modes are similar 
for both classes of models. For models in which the net magnetic 
flux vanishes initially we find that M^™ is roughly constant for 
sufficiently strong initial fields, the stress being slightly larger 
for sinusoidal (^l^p) than for step-like initial fields (^~tep)- 
models develop a more complex field morphology with more 
intense current sheets and more potential sites for reconnection 
than uniform field models. Thus, the growth rates of the resistive 
instabilities are comparable to those of the MRI for much weaker 
fields than in the uniform-field models. MRI termination also 
occurs at smaller values of Mj^™. 

4.3. Axisymmetric models with entropy gradients 

We also simulated some axisymmetric models imposing an ad- 
ditional entropy gradient. In this case, the instability criterion 
is more complicated and various instability regimes exist (see 
Sect.|2]): (magnetic) shear instability, convection, and magneto- 
buoyant instability. Otherwise unstable modes may be stabilized 
by a stable thermal stratification or by fast (not necessarily differ- 
ential) rotation. Tables lBTI (models with a positive entropy gra- 
dient) and |B.2| (models with a negative entropy gradient) provide 
a list of the simulated models. Note that all models discussed in 
this subsection have a uniform initial magnetic field. 



4.3.1. Positive entropy gradients 

We first discuss differentially rotating models having a stabiliz- 
ing entropy gradient comparing models with and without mag- 
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Fig. 14. Maxwell stress M^™ at MRI termination as a function 
of the initial magnetic field strength h'^ for models with non- 
uniform initial fields. Models with and without velocity damp- 
ing are shown by asterisks and diamonds, respectively. The blue, 
green and red symbols denote models where the z-component of 
the initial magnetic field is given by ^l^p (see Eq. (l44ll). h\^^^ (see 
Eg. (|46]l), and h\., (see Eq. (|45]l), respectively. The light gray 



lines illustrate power laws oc ^^7~^ 
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netic field. The non-magnetic models are stable due to the posi- 
tive entropy gradient, i.e. initial perturbations do not grow. 

The models with a magnetic field belong to the MSI regime 
(cf. Fig.|2]i, their MRI growth rates being reduced compared to 
models with no entropy gradient. We simulated models with dif- 
ferent entropy gradients (d^S = 0.02, 0.04, 0.08 km"') and dif- 
ferent adiabatic index of the equation of state (Fb - 1.31,5/3). 
Generally, we find a good agreement between the analytic pre- 
dictions and the numerical results. For dj^S < 0.08 the mod- 
els are unstable belonging to the MSI regime, whereas an en- 
tropy gradient of dj^S - 0.08 suffices to stabilize the model. The 
growth rates agree well with the analytic ones, and the numer- 
ical models show the typical dependence of the growth rate of 
the MRI on the initial magnetic field strength. However, there 
exists one interesting difference: cr increases from small values 
for weak initial fields for which /Imri is under resolved and con- 
verges to the correct growth rate for strong fields for which /Imri 
exceeds the grid resolution significantly. Unlike for models with- 
out entropy gradient, the growth rate becomes largest for mag- 
netic fields for which the MRI wavelength is similar to the grid 
resolution, and at these field strengths the numerical growth rates 
can exceed the theoretical ones. 

Dynamically the models behave similarly as models without 
an entropy gradient. Channel flows develop during the growth 
phase of the MRI, their width being set by the MRI wavelength. 
MRI termination occurs due to the growth of tearing-mode-like 
resistive instabilities. When velocity damping is applied, the 
maximum Maxwell stress at MRI termination, M"^™, is deter- 
mined by the box size. Comparing models with a positive en- 
tropy gradient and with no entropy gradient we find a common 
linear relation between M^™ and ctmri, indicating a common 
reason for MRI termination (see Fig.fTsTl. 

According to Eq. ( l28T l, the channels are wider in models with 
larger entropy gradients. As wider channels are less prone to re- 



M. Obergaulinger et al.: MRI in core collapse supernovae 



15 





4 


; 














3 












- 
















o 
















o 


2 










- -o- - 




s 














1 





































O 









0.20 0.40 0.60 0.80 1 1.20 

Fig. 15. Maxwell stress at MRI termination, M^™, as a function 
of the MRI growth rate, o-f^jf^i, for models with zero (diamonds) 
and positive (plus signs) entropy gradients. 



sistive instabilities, they can support stronger fields before be- 
ing disrupted (see discussion above), i.e. M^™ is larger in mod- 
els with larger entropy gradients. The MRI growth rate, on the 
other hand, is smaller for larger entropy gradients (see Eq. ( |29] l) 
implying that Mj^™ decreases with increasing entropy gradient. 
Both effects taken together suggest a weak anti-correlation of 
Mj^™ with the size of the (positive) entropy gradient. An anti- 
correlation is also suggested by our numerical results, although 
more models are needed to confirm it. It is unclear, for exam- 
ple, whether the growth rates of resistive instabilities derived in 
Sect. 14.2.21 also hold for stably stratified media, and whether the 
boundary conditions have an influence in models with large en- 
tropy gradients. Small perturbations of the quasi-periodic (be- 
cause of global gradients) radial entropy distribution may leave 
their imprint on MRI termination by enforcing a preferred length 
scale, thus clouding effects due to a variation of dj^S . 

We have also simulated a few of the models using an ideal- 
gas equation of state, P = (F - l)e, instead of the hybrid EOS 
finding, however, no effect on the evolution of the models. 

4.3.2. Negative entropy gradients 

Convection can develop in models having a negative entropy gra- 
dient, but it can be suppressed by rapid rigid or differential rota- 
tion. If a magnetic field is added to a convectively unstable sys- 
tem which cannot be stabilized by rotation, the system is located 
in the "convection" regime of Fig.|2] In a situation where rota- 
tion suffices to suppress the convective instability, the addition 
of a weak magnetic field puts the system into the "MBI" regime, 
and a magneto-buoyant instability similar to the standard (i.e., 
d^S = 0) MRI develops. When d^S < the equivalent of the 
standard MRI corresponds to the "MSI" regime in Fig.|2l 

Before discussing our results (see Tab. |B.2| for a list of the 
simulated models), we need to comment on the boundary con- 
ditions. Allowing for radial transport of energy shearing-disc 
boundaries can, in principle, lead to a transport of entropy across 
the pseudo-periodic boundaries, thus modifying the initial en- 
tropy profile. By comparing shearing-disc models and models 
with reflecting boundary conditions, we verified that none of 



these boundary conditions suppresses the growth of the MRI, 
and that both give similar results. 

Let us first consider a non-magnetic model which rotates 
rigidly with an angular velocity Qo = 1000 s"' and has an en- 
tropy profile given by So - 0.2 and dj^S - -0.075 km"' (see 
Eq. dSTll). With ~ -14 and H^/nl = 0, the model be- 

longs to the convective regime. Buoyant modes are unstable 
and grow at a theoretical rate o-± = 3.4 ms"'. Simulated on a 
grid with 2.5 m spatial resolution, the model is unstable with 
a numerical growth rate cr - 2.7 ms"', and convection sets in 
quickly. The flow is dominated by a few (one or two) fairly cir- 
cular convective rolls. Due to the transport of entropy and an- 
gular momentum by the overturning fluid, the model develops 
complementary entropy and rotation profiles characterized by 
"cold" (i.e., low-entropy), rapidly rotating matter in down-flows 
and "hot" (i.e., high-entropy), slowly rotating matter in up-flows. 
The redistribution of angular momentum and entropy leads to an 
average (with respect to the z-coordinate) rotational profile of 
the form Q. oc cr"^, i.e., constant specific angular momentum 
(see Fig.fTSTl. and a flat entropy profile. 

For a faster rotation rate of f2o - 1500 s"', corresponding to 
N^/Q.^ ~ -5.7, i.e., still in the convective regime, the evolution 
is similar except for a reduced growth rate (cr a: 1.6 ms"') due to 
rotational stabilization. The model develops differential rotation 
with the same cr-dependence as in the case of slower rotation. If 
we increase the rotation rate to Qo = 1900 s"' (A^^/Qq x -3.2) 
buoyant modes are stabilized by rotation. 

The above results also hold if the initial model is rotating 
differentially. In particular, convection (i.e., the negative initial 
entropy gradient) gives also rise to a rotation law of the form 
Q oc zir"^, and a flat entropy profile. 
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Fig. 16. Q, averaged over z, as a function of m for non-magnetic 
convective models at f = 5 1 .9 ms (green dashed), with an initial 
magnetic field b'^ - 10'^ G at f = 51.9 ms (red dash-dotted), and 
an initial field - lO'-' G (blue solid lines), respectively. For the 
latter model the different symbols indicate different epochs: f = 
11.5 ms (plus signs), t - 20.8 ms (diamonds), and f - 41.7 ms 
(squares). The dotted black lines show the initial rotation law 
(Qo = 1500 s"', and a power-law profile Q oc tir"^. 

Next, we add a magnetic field of - lO'^G to a convective 
model, i.e., to a model with a negative entropy gradient rotating 
too slowly (in our case, for rigid rotation, with Qo ^ 1500s"') 



16 



M. Obergaulinger et al.: MRI in core collapse supernovae 



28 



^ 27 
> 



o 

24 



23 



22 















- 


















f 

I 








- -{ 

1 
( 
(' 








T 
t 








1 
f 

it 









10 20 30 

t [ms] 



Fig. 17. Evolution of the mean magnetic energy density e™^^ 
(solid black line), the mean energy densities corresponding to 
the m (dotted red), (dashed brown), and z (dash-dotted green) 
component of the magnetic field, and the absolute value of the 
mean Maxwell stress component Mnj,^ (dashed blue line) for the 
model whose rotational profiles at different times are shown in 
Fig-HSl i.e., a rigidly rotating model (Qq = 1500 s"') with an 
initial magnetic field of Z?;, = 10'^ G. 



for convection to be stabilised by rotation. The temporal evo- 
lution of the magnetic energy and the Maxwell stress of this 
model is shown in Fig. [17] while Fig. [18] gives the spatial dis- 
tribution of the entropy of the model at two different times. 
Initial perturbations are amplified rapidly, but saturation sets in 
within 6 milliseconds the growth rate being slightly higher than 
in the non-magnetic models. With N^/Q.^ » -5.7, the model 
is dominated by buoyant modes, but there still exists some in- 
fluence of the Alfven modes. In particular, although infinitely 
long modes grow rapidly, the fastest growing modes are Alfven 
modes of finite wavenumber which depends on b'^. For suffi- 
ciently strong initial fields (or sufficiently fine resolution), these 
modes are numerically resolved, and have a growth rate exceed- 
ing that of the corresponding non-magnetic model. The magnetic 
field strength increases exponentially as the instability devel- 
ops, and at the onset of saturation large convective rolls develop 
(left panel of Fig. [18]). In the saturation phase (right panel of 
Fig.[T8]l the flow geometry differs considerably from that of the 
corresponding non-magnetic model. It consists of down-drafts 
of cold material and up-flows of hot gas forming small-scale 
structures rather than large circular convective rolls. Like in the 
non-magnetic model, cold and hot regions correspond to regions 
of low and high angular velocity, respectively. Differential rota- 
tion with constant specific angular momentum, Q oc m^^, de- 
velops due to hydrodynamic transport of angular momentum in 
convective overturns. The magnetic energy related to the radial 
component of the magnetic field and the Maxwell stress compo- 
nent Mjjr^ remain high during saturation, i.e. angular momentum 
transport converts the m^^ rotation law prevailing at early epochs 
into nearly rigid rotation (Fig. [16]). 

If for the same initially rigidly rotating model the initial mag- 
netic field is too weak to resolve /Imri (we simulated two mod- 
els with = 10"^ and 10'^, respectively) convection develops. 
The growth rates are similar to those of non-magnetic models. 
The weakest initial field, fejj - 10^*0, has no impact at all, and 



the evolution of the model with an initial field of 10 G dilTers 
only slightly from that of the non-magnetic model. The magnetic 
field increases exponentially as the instability grows. Large per- 
sistent convective rolls form, and differential rotation develops. 
After an initial exponential growth the mean magnetic energy re- 
mains large, but the contribution of the radial component of the 
magnetic field and the mean Maxwell stress component M^^ de- 
crease almost by four and two orders of magnitude, respectively. 
Consequently, no significant angular momentum transport oc- 
curs due to magnetic stresses, and similar to the non-magnetized 
model a Q oc m^^ rotation law develops. Fig.[T6b. 

We now consider the MBI regime (see Fig.|2| and discuss 
models rotating initially rigidly with Qq = 1900 s ' and hav- 
ing an initial entropy gradient d^jS - -0.1 km ' (see Tab. |B.2t , 
which implies A^^/Qg ~ -3.6. Without magnetic field, the insta- 
bility of the buoyant modes is suppressed by the fast rotation. 
However, if a weak magnetic field is added, an instability of the 
MBI type develops, i.e., the Alfven modes become unstable. The 
numerical growth rates show a similar dependence on the mag- 
netic field strength as in case of the standard (i.e., d^S = 0) 
MRI, because Amri is resolved. The instability grows rapidly 
(cr ~ 1.4ms"', similar to the theoretical value crih x 1.7 ms"'). 
During the growth phase channel modes appear, which lead to a 
transport of both angular momentum and entropy. After an expo- 
nential initial growth and some decrease after MRI termination 
the mean magnetic energies contained in the total magnetic field 
and all three field components remain large (corresponding to 
field strengths of ~ lO'^'G), but the mean Maxwell stress com- 
ponent Mjjj^ drops to zero within ten milliseconds oscillating af- 
terward with decreasing amplitude between positive and nega- 
tive values. Hence, large-scale angular momentum transport is 
limited. At the end of the simulation, the model shows consid- 
erable variations in Q., but there is no clear indication of a mean 
dilTerential rotation of the form Q. = Q(tiT). The entropy profile 
in the saturated state is almost flat. 

Finally, we summarize a few common features of the mod- 
els having a negative entropy gradient (see Tab. |B.2| l. All models 
develop instabilities in accordance with the flow regime to be 
expected from their model parameters. The growth rates of all 
models are, within the uncertainties, similar to the theoretical 
predictions. As for the dynamics, we have to distinguish mod- 
els in the convective regime from those in the mixed and MBI 
regimes. The former class of models shows convective mush- 
rooms and large-scale overturns with only little influence of the 
magnetic field, whereas the last class is dominated by channel 
flows. Consequently, angular momentum transport by hydrody- 
namic flow leads to a rotational profile Q - vj"^ for models in 
the convective regime, while models in the mixed-type regime 
tends towards rigid rotation the angular momentum transport 
being dominated by magnetic fields. Termination of the insta- 
bility growth occurs for models both in the MBI and mixed-type 
regime analogously to that of models in the MSI regime without 
entropy gradient, i.e., by reconnection in resistive instabilities al- 
tering the topology of the channel flows. Consequently, we find 
similar dependencies on the initial field strength, the grid reso- 
lution, and the type of boundary conditions. The instability in 
convective models, on the other hand, saturates when the initial 
entropy gradient is removed by vigorous entropy transport due 
to overturning fluid motions. 

4.4. Three-dimensional models 

The results of the axisymmetric simulations discussed in the pre- 
vious section demonstrate the possibility of MRI-driven field 
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Fig. 18. The entropy distribution (color coded), the poloidal velocity field (arrows), and the magnetic field lines of the model 
displayed in Fig.[T7]at t = 7.3 ms (left panel) and f - 20.8 ms (right panel), respectively, i.e., at the begin of the saturation phase 
and during saturation. 



amplification in core collapse supernovae, and provide some 
insight into the evolution of MRI unstable layers in the core. 
However, to address the MRI problem in full generality, one has 
to consider three-dimensional models, because the assumption 
of axisymmetry implies severe restrictions for the dynamics of 
the magnetic and kinetic fields. The most important limitations 
are that, in axisymmetry, a toroidal field cannot be converted into 
a poloidal one, and that the disruption of the c hannel flows re- 
quire s non-axisymmetric parasitic instabilities (iGoodman & Xul 
11994 . 

As 3D simulations are computationally much more expen- 
sive than 2D ones, we could not perform a comprehensive study, 
but had to focus on a few selected models. We simulated mod- 
els with different field geometries and varied the initial field 
strength, the entropy profile, and the grid size (see Tab. |B.3| and 
Tab.ltHi. 

4.4.1. Uniform initial magnetic fields, no entropy gradient 

We first discuss models which have a uniform initial magnetic 
field bf^ in z-direction, no entropy gradient, and rotate differen- 
tially with Qo = 1900 s-i and an = -1.25 (see Eq. ((33). If the 
MRI wavelength is well resolved (e.g., for models with initial 
field strengths of 2x lO'"* G and 4x lO'"* G simulated at a grid res- 
olution of Snr - 20 m), the growth rate is high and independent 
of . Under-resolved models (e.g., models with bf^ = 10'^ G 
simulated at the same resolution) exhibit a slower growth of the 
MRI. From the growth rates of the MRI, we infer that in 3D the 
same resolution criterion applies as in the case of axisymmetry. 

During early epochs the evolution is similar to that of the cor- 
responding axisymmetric models: a number of radially aligned 
channels appear. Strong differential rotation causes significant 
wind-up of flow features leading to structures elongated in 0- 
direction, i.e., there exists only a modest variation of the MHD 
variables with azimuthal angle at this stage. Sheet-like struc- 
tures dominate the field geometry. The rotational profile begins 
to show distortions due to the transport of angular momentum 
by Maxwell stresses (see left panel of Fig.[T9]l. At later epochs 
the flow in 3D is more complex than in axisymmetry. Although 
coherent structures, i.e., flux sheets, are still present, their geom- 



etry is more tangled and twisted, and less isotropic than earlier 
in the evolution (see middle panel of Fig.[T9]l. 

An evolution from coherent channel flows to a more tur- 
bulent state is characteristic for all three-dimensional models 
with a uniform initial magnetic field. However, as pointed out 
by Sano & In utsukal (1200 II) . channel flows can develop again 
from the turbulent state. Consequently, the magnetic field can 
continue growing, and the angular momentum transport will 
be enhanced strongly. In the most extreme cases, the evolu- 
tion is similar to that of a corresponding axisymmetric model. 
This is exactly what we observe for some models at late times, 
t > 30 ms (see right panel of Fig.[T9]l. when a dominant chan- 
nel flow forms. These model enter again a state of exponential 
growth, and a large part of the angular momentum is extracted by 
Maxwell stresses. The field strengths reach several 10'^ G, peak- 
ing at lO'^ G, and the mean Maxwell stress component M^^ ex- 
ceeds 10^° erg cm"-' (see middle panel of Fig.l20]l. and compare 
with the corresponding axisymmetric model in the left panel). 
Despite a qualitative similarity between the evolution of the 3D 
and axisymmetric models, we note that the secondary exponen- 
tial growth is slower in three dimensions. 

The emergence of a large-scale structure of the magnetic 
field from a turbulent state can be seen in Fig.|2T] comparing 
the field structure at f = 26.8ms and t - 42.5ms, respectively. 
At the earlier time (left panel), we find a small-scale field dom- 
inated by slender flux tubes. Field lines of different polarity (in- 
dicated by different colors) are lying close to each other. After 
the development of the channel flow (right panel), the field is 
dominated by a large-scale pattern. A smooth surface perme- 
ating the box at nearly constant z-coordinate separates in two 
large regions field lines of different polarity from each other. In 
each of the two regions, we find one broad flux sheet where most 
of the magnetic energy is concentrated. The separation layer is 
filled by gas rotating nearly uniformly at a ow angular veloc- 
ity (Q. ~ 1500s"'). The surrounding gas rotates uniformly as 
well, but at a much higher velocity (Q. ~ 1800s"'). The two 
flux sheets form a thin transition region between both rotational 
states. Thus, the dynamics is similar to that of the corresponding 
axisymmetric model. 
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Fig. 19. Structure of a 3D model with - 2x 10^^ G computed on a grid of 1 km^ at a resolution of 20 m at f = 16.2 ms (left), 
t - 26.8 ms (middle), and t = 42.5 ms (right), respectively. Shown is the volume rendered magnetic field strength (blue to green), 
and a red-orange iso-Q surface corresponding to O = 1820 s ' (left and middle) and Q. = 1680 s ' (right), respectively. The red, 
green, and blue axes point into m, 0, and z direction, respectively. Channel flows can be identified as green sheet-like structures. 




Fig. 20. Evolution of the mean magnetic energy density e™^^ (solid black line), the mean energy densities corresponding to the 
m (dotted red), (dashed brown), and z (dash-dotted green) component of the magnetic field, and the absolute value of the mean 
Maxwell stress component M^^^ (dashed blue line) for models with an initially uniform magnetic field b^^ = 2x lO'-' G in z-direction, 

and a rotation law given by Qo = 1900s ' and an - -1.25. The panels show a 2D model computed in a box of LnjXi; - Ikmxlkm 
(left), a 3D model computed in a box of Lj^ x x L- - 1km x 1km x 1km (middle), and a 3D model computed in a box of 
Lj;j X L^x L, - 1km x 2km x 1km (right), respectively. The grid resolution is 20m in all three cases. 



Because our boundary conditions allow for a loss of angular 
momentum, and thus for the total disruption of the differential 
rotation profile by transport through the radial boundaries, this 
stage represents the end of the evolution, just as it did in ax- 
isymmetry: the instability has used up its free-energy reservoir 
Hence, the later evolution consists only of violent oscillations. 

Only a subset of our models show a prominent re-appearance 
of single channel flows, and most of them do not exhibit a sec- 
ondary exponential growth phase. Instead, the mean magnetic 
energy and the Maxwell stress remain roughly constant dur- 
ing saturation, albeit fluctuating strongly (see the right panel 
of Fig.|20t. Angular momentum transport is less efficient for 
these models, and their initial rotation profiles remain nearly 
unchanged. A turbulent flow and magnetic field persist during 
saturation, and coherent, channel-like structures develop tran- 
siently. The structure of the magnetic field of a model with 

= 4 x lO'^G computed on a grid of 50 x 100 x 50 zones is 
displayed at two different epochs in Fig.|22l At f = 21.5 ms (left 
panel) one recognizes a turbulent state, while large-scale patterns 
(right panel; yellow structures) dominate the flow at f = 37.2ms 
when the magnetic field strength is largest, and the Maxwell 



stress is strongest. Unlike in the model discussed above, the co- 
herent flow is unstable and becomes turbulent within a few mil- 
liseconds, and the absolute value of the Maxwell stress \m^^\ 
decreases. 



Channel modes and parasitic instabilities: The appearance 
and stability of single large-scale flows that lead to a secondary 
exponential growth phase and eventually to the disruption of the 
rotation profile depend on the geometry of the simulated domain 
as well as on the ratio of the grid resolution and the fastest grow- 
ing mode. 

Models, which are computed in a box of Ikm^ with a res- 
olution of 20m and where velocity damping is applied, develop 
secondary stable channels, if the initial magnetic field is stronger 
than 2 x lO'^G. The MRI growth rates found for these models 
(cr = {0.76, 1.03, 1.10) ms-i for bo = {1,2,4} X lO'^G, respec- 
tively) indicate that the grid resolution is sufficiently fine to re- 
solve the fastest growing MRI mode for the two most strongly 
magnetized models. However, it is too coarse for the model with 
the weakest initial field, because the theoretical growth rate for 
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Fig. 21. Same as middle and right panel of Fig. [19] but showing besides the volume rendered magnetic field strength (blue to green) 
also the magnetic field lines, which are obtained by starting the integration of the magnetic field at two surfaces of constant 0- 
coordinate (i.e., orthogonal to the green axis) at the left and right hand side of the domain. The field lines originating from the 
left and right surface are plotted in red and yellow, respectively. The right panel shows, in addition, the isosurface - (i.e. the 
magneto-pause). 




Fig. 22. Volume rendered magnetic field strength of a model with b^^ - 4 x lO'^ G computed in a box of 1 X 2 X 1 km^ with a 
resolution of 20m at f = 21 .5ms (left) and t - 37.2ms (right), respectively. The coordinate directions are indicated as in Fig. [19] 



the fastest growing MRI mode is ctmri = 1.14 ms ' for these 
models (see Sect. 14. 2. Ti l. 

To investigate the stability properties of large-scale channel 
modes as a function of the box geometry, we simulated mod- 
els with an initially uniform magnetic field using boxes of dif- 
ferent size and shape. The models were rotating according to 
Qo = 1900 s"' and an - -1.25, and their initial magnetic field 
was fog = 4 X lO'^G when applying velocity damping bound- 
aries, and b^-lx lO'^G otherwise. We varied both the ratio be- 
tween the radial and vertical, L^/ic and the radial and toroidal 
size , L^/ Lgj, size of the box. The grid resolution was 20 m (see 
Tab.|B3]). Plotting the stress ratios M'^flM'^J (Fig.|23l damp- 
ing boundaries) and {M^^)IM^^ (Fig.|24l non-damping bound- 
aries) as a function of the aspect ratio of the computational box, 
provides some indication of the range of M^t^ values prevail- 
ing during the post-growth phase. The ratios allow one to dis- 
tinguish models with a strong variability due to the dominant 
re-appearance of channel modes from those models exhibiting a 



smooth evolution without dominant large-scale coherent struc- 
tures. 

We find that models with a radial aspect ratio Lj^lL- - 1 
and a toroidal aspect ratio L^/L^ > 2 are unstable against par- 
asitic instabilities, independent of the grid resolution in toroidal 
direction. Turbulence develops and leads to a flow structure as 
shown in Fig. [22] Models having the same radial aspect ratio, 
but a smaller toroidal one are stable and evolve similarly as ax- 
isymmetric models, i.e., parasitic instabilities do not grow and 
a dominant large-scale channel flow develops, which gives rise 
to a morphology of the type presented in Fig. [22] These find- 
ings do not depend on how the growth of the MRI ends, i.e., 
whether velocity damping is applied and reconnection between 
adjacent channels occurs inside the box, or whether no damp- 
ing is imposed and reconnection occurs near the surface of the 
computational box. 

These results can be understood from th e analysis of par- 
asitic instabilities by iGoodman & Xul (il994t) . who argued that 
three-dimensional flows are unstable against parasitic instabili- 
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ties, but these instabilities can be suppressed by the geometry of 
the computational box. According to their analysis, the growth 
rate of the parasitic instabilities is highest for modes with half 
the wave number of the unstable MRI modes they are feeding 
off. Hence, if a channel flow forms at late times with a wave- 
length equal to the box size in z-direction, L^, unstable parasitic 
modes must have a toroidal wavelength ~ 2L~ to grow rapidly. 
Thus results in the criterion for the channel flow instability we 
have found in our simulations. 

In a ccordance with simulations presented recently by 
iBodo et al. (2008), we find that models with a radial aspect ra- 
tio Lcj/L, < 1 experience a second exponential growth phase as 
described in Sect.|42] (note the large ratios of M™" and Mj^™ 
for the corresponding models in Fig.l23ll. whereas a larger ra- 
dial aspect ratio appears to favor a less violent post-growth phase 
where coher ent channel modes c an appear but are disrupted after 
a short time. iBodo et al.1 (l2008l) obtained this result for simula- 
tions performed with a toroidal aspect ratio L^jL. - 4. 

We confirm a similar dependence of the dynamics on the ra- 
dial aspect ratio also in axisymmetry and for three-dimensional 
boxes with smaller L^/L.. In this case, parasitic instabilities are 
unable to disrupt the channel modes. Consequently, the MRI 
experiences a second exponential growth phase dominated by 
just one (two in a few cases) large channel mode of width a 
which is determined by the size of the computational box in 
z-direction. The maximum Maxwell stress that can be reached 
is limited by the onset of resistive instabilities. The dependence 
on the channel width (Eq. (lA.llb ) explains why the maximum 
Maxwell stress varies with L.: larger boxes allow for wider chan- 
nels for which the resistive instabilities grow slower, thus requir- 
ing higher Alfven velocities for a growth rate comparable to the 
one of the MRI. Hence, the MRI reaches stronger fields for larger 
(in z direction) boxes. 

Despite the differences in the MRI termination process, the 
behavior of models with and without velocity damping is quite 
similar, because the velocity damping does not affect the sec- 
ond generation of vigorous channel flows significantly. Thus, the 
breakup of these channels and the values of the corresponding 
maximum Maxwell stress do not depend strongly on the choice 
of the boundary condition. On the other hand, MRI termination 
(the termination of the initial exponential growth of the MRI) 
does depend on whether velocity damping is applied or not, and 
thus the ratio M^^/M^™, too. 

For boxes having a large toroidal aspect ratio, L^/L, > 2, we 
observe a quiet evolution during the non-linear saturation phase 
of the MRI when varying the radial aspect ratio, L^/L.. In par- 
ticular, the fluctuations of M^jj^ are small after MRI termination 
for models where both aspect ratios are large. The models in the 
upper right corner of Fig. |23] have values of {M^^) jM^^^ close 
to unity. 

We may try to infer some consequences from these results 
for the MRI in supernova cores. Close to the core's equator the 
region, where the MRI develops, can have to a small radial size, 
L^, ranging from a few to a few hundred kilometers determined 
by the gradients of O and S in the core. The vertical extent of 
the unstable region, L., can be expected to be of similar size. 
The azimuthal extent of the MRI unstable region, L^, will be 
significantly larger, leading to a non-violent evolution of the sat- 
urated state of the MRI. The geometry is different close to the 
pole. However, we cannot apply our results there without modi- 
fications as we have considered only cases where the gradients of 
f2 and of all thermodynamic quantities are aligned - a situation 
which does not apply near the poles. 
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Fig. 23. The left panel shows the ratio of the maximum Maxwell 
stress per unit volume, M^^, and its value at MRI termina- 
tion, M|^™, as a function of the toroidal and radial aspect ra- 
tios, L^/L, and L^/L. for the models listed in Tab. IB. 3 1 The right 
panel shows the ratio of Mj;j^ averaged over the saturation phase 
and the value at MRI termination. Each model is represented by 
a symbol its color reflecting its maximum Maxwell stress. All 
models are computed imposing velocity damping at the radial 
grid boundaries. 
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Fig. 24. Same as Fig.|23] but for models where no velocity damp- 
ing is applied at the radial grid boundaries. 



Effects of resolution and initial magnetic fields in 3D vs 2D: 

After having discussed how the aspect ratios of the simulation 
box determine the Maxwell stress at MRI termination and during 
the subsequent saturation phase, we will now compare whether 
the behavior of 3D models differs from that of 2D models. The 
models discussed in this paragraph are listed in Tab. |B.3| and |B.4| 

First, we note that in 3D, as in axisymmetry, the growth rate 
of the instability is not affected by the choice of the grid provided 
the fastest growing mode is resolved. 

Fig.|25]demonstrates that in 3D the dependence of M^™ on 
the initial magnetic field strength is well described by the same 
power law as in axisymmetry: models without damping of the 
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Fig. 25. Volume-averaged Maxwell stress component M^™ at 
MRI termination as a function of the initial magnetic field 
strength, bo, for axisymmetric models with uniform initial mag- 
netic field in z-direction, a rotational profile f2 - 1900 s ' cr 
and a vanishing entropy gradient for a set of 3D models com- 
puted with (asterisks) and without (diamonds) velocity damping 
at the radial boundaries. Models computed in a box with a radial 
size of 1 km and 0.5 km are shown with green and blue symbols, 
respectively. The colored bands are the same as in Fig.[8] 

radial velocity line up along a band oc (Zjjj)'^''^, and models with 
velocity damping are characterized by a roughly constant value 
of Mj^™, which depends on the size of the radial box (compare 
with Fig. [8]). This agreement is to be expected as the growth and 
the resistive disruption of channel flows are essentially axisym- 
metric processes, which are, thus, not significantly modified by 
three-dimensional effects. 

After MRI termination the evolution of M^^ depends on the 
aspect ratio of the computational box (see discussion above). 
When averaging the fluctuating Maxwell stresses over the sat- 
uration phase, we find values for (Mcj^) which differ consider- 
ably from those of Mj^™. Lacking a thorough understanding of 
the instabilities involved in the MRI saturation process, and hav- 
ing only a imited set of 3D models at hand, one is not yet in a 
position to formulate a better description of the dependence of 
the evolution after MRI termination on the aspect ratio of the 
box, and to provide a unified description of MRI saturation am- 
plitudes. 

4.4.2. Uniform field, entropy gradients 

Mixed regime: Let us first consider models from the mixed 
regime having an initial rotation law given by Q.o - 1900 s"' 
and an - -1-25, and an entropy distribution given by So - 0.2 
and djjjS = -0.038 (i.e, with a negative entropy gradient). In ax- 
isymmetry the MRI grows in these models with a rate of ctmri ~ 
1.7ms"', i.e., close to the theoretical value of x 2.0ms"'. The 
3D models show the same growth rate provided the spatial grid 
resolution is sufliciently high. 

The long-term evolution (i.e., many rotational periods into 
the non-linear phase) of the models depends strongly on the 
choice of the radial boundary conditions. If the entropy at the 
inner and outer boundary is allowed to change (i.e., using reflec- 
tive boundaries), a flat entropy profile develops after a short time. 
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Fig. 26. Evolution of the mean magnetic energy density e™° 
(solid black line), the mean energy densities corresponding to 
the Tu (dotted red), (dashed brown), and z (dash-dotted green) 
component of the magnetic field, and the absolute value of the 
mean Maxwell stress component M^.^j (dashed blue line) for a 
3D model belonging to the mixed regime. The model rotates dif- 
ferentially with Qo = 1900 s"' and an = -1.25. The initial en- 
tropy gradient is d-u,S - -0.038 km"', and the initial magnetic 
field strength is = 2 x lO'^G. The model was simulated in 

1x2x1 km""' and on a grid of 



a box of size L-^j x L^x L, 
50 X 100 X 50 zones. 



To reduce the influence of boundary effects, one could employ a 
technique widely used in simulations of convective layers: add 
a cooling layer on top of and an overshoot layer below the con- 
vection zone. However, exploring this approach was beyond the 
scope of the present work. 

For models having a negative entropy gradient the growth 
of the MRI is not influenced by 3D effects, if the fastest grow- 
ing mode is resolved. Thus, their behavior is similar to that of 
models with no entropy gradient. The Maxwell stress at MRI 
termination also does not differ significantly from that of the 
corresponding axisymmetric models, and due to the boundary 
conditions applied in the models (velocity damping) its value is 
set by recombination of field lines close to the inner and outer 
radial boundary. 

Contrary to axisymmetric models, the saturated state of the 
3D models does not show any sign of a late exponential growth 
phase characterized by the re-appearance of channel modes, and 
the saturated MRI stresses are smaller in magnitude than M^™, 
i.e., the maximum Maxwell stress is reached at MRI termination 
(see Fig.|26land compare with Fig.|4]i. The evolution of the av- 
erage radial entropy profile profile, computed as the average of 
S (m, (p, z) at constant m, is shown in Fig.|22] Until the saturation 
of the instability (at t » 11ms), the initial linear profile S(-U7) 
is basically unchanged. However, afterward the entropy profile 
flattens, S becoming nearly constant for 15.2km < tir < 15.8km. 
Close to both radial boundaries, the entropy profile develops ex- 
trema, which are most likely an artifact of our boundary condi- 
tions. The flat entropy profile is stable and does not vary strongly 
with time. The Q. profile flattens after the initial growth phase, 
too. The velocity field in the saturated state is dominated by a 
rich small-scale structure, while the magnetic field is organized 
in a multitude of flux tubes. 
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Fig. 27. Average radial entropy profile as a function of time for 
the model shown in Fig.l26l 



MBI regime: Next we consider a few models that, in axisym- 
metry, belong to the MBI regime. Initially, the models rotate 
rigidly with Q = 1900 s"' and possess an entropy gradient 
d-cjS - -0.10km"' (i.e., C - -3.6). We computed models for 
bl = 10 G, 10'°G, and 2 x lO'^G, respectively. All models ai-e 
simulated in a box of size L^jXi^xL, = lx2xxl km"^ and 
on a a grid of 50 x 100 x 50 zones. 

Contrary to their axisymmetric counterparts (see Sect. l4.2l i. 
these models develop convective modes e ven when no m ag- 
netic field is present. As described, e.g., bv iTassouH (Il978h . ro- 
tation can stabilize axisymmetric modes in a convectively sta- 
ble environment, but non-axisymmetric modes can nevertheless 
grow in that situation. The model with the weakest initial field 
(Z7g = lOG) shows a growth of non-axisymmetric MBI modes, 
but these modes cannot be resolved due to the extremely weak 
initial field. Thus, the model behaves essentially similar to an un- 
magnetized one, but can serve as a reference model for initially 
more strongly magnetized models where we can resolve /(mri- 
The growing convective modes eventually extend over the en- 
tire domain in radial and z-direction, while having a small wave- 
length in (f) direction (see Fig.|28] upper panel). The exponen- 
tial growth of the convective instability saturates at f =s 7ms. 
During this growth phase the mean magnetic energy increases at 
the same rate as does the kinetic energy. After MRI termination 
the structure of the model is characterized by two large, roughly 
cubic convective cells with a size of about 1 km"* instead of a 
multitude of elongated structures (see (Fig.|28] lower panel), and 
an essentially flat entropy profile. The magnetic field is subject 
to kinematic amplification at a smaller growth rate than before 
MRI termination due to stretching in the convective vortices. At 
much later epochs the typical size of structures in the velocity 
field decreases again, leading to more turbulent fields. The mag- 
netic field, which is too weak to affect the dynamics, is passively 
advected with the flow. 

For the model having the strongest initial magnetic field 
b^Q - 2 X lO'-' G we can resolve /Imri- The axisymmetric ver- 
sion of this model showed a MBI growth a rate close to the 
theoretical one (ctmri ~ 1.7 ms"'). For the 3D model we find 
Cconv ~ 2.6 ms"'), i.e., its evolution is dominated by convection 
(although we are able to resolve the MBI), and the MBI growth 
rate is similar to that of a weakly magnetized model (see previ- 




Fig.28. Flow structure of a MBI model with b^ - lOG at 
t - 4.86ms (upper panel), and t - 14.4ms (lower panel), re- 
spectively. The solid black lines are stream lines of the velocity 
field computed in a frame co-rotating with the mean angular ve- 
locity, and regions of positive and negative radial velocity are 
colored in red and green, respectively. The colored arrows point 
into the same coordinate directions as in Fig.[T9l 



ous paragraph). MBI growth is mediated by non-axisymmetric 
modes having the same elongated geometry as those in the es- 
sentially unmagnetized model. After saturation, a few large vor- 
tices of approximately cubic shape form, which later decay into 
small-scale structures again. An intermediate stage of this de- 
cay process is displayed in Fig.|29] when one large vortex is still 
present in the right half of the box, while its left half is dom- 
inated by spatially less coherent fields. At even later times the 
vortex disappears and the structure of the whole model is simi- 
lar to that shown in the left half of the box. The mean magnetic 
energy and Maxwell stress are small compared to typical MSI 
or mixed models. Compared to a differentially rotating model 
(Q.Q - 1900 s"', an - -1.25) with a vanishing entropy gradient, 
the maximum magnetic fields are reduced by a factor of < 2, and 
the mean magnetic energies and Maxwell stresses by a factor of 
~ 10, but we still find a slow growth at the end of the simulation. 

Finally, we add a few comments on a model where we cannot 
resolve Amri (b^ = 10'"G), but where the magnetic field satu- 
rates within ~ 60ms after the onset of convection. The model 
evolves similarly to the essentially unmagnetized one, but at 
t X 60ms the energy of the kinematically amplified magnetic 
field becomes almost as high as the convective kinetic energy. 
The amplification process ceases, and the magnetic energy lev- 
els off'. Close to the end of the simulation convective transport 
gives rise to a rotation law Q x const.. 
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Fig. 29. Structure of a rigidly rotating model with /7~ = 2 x lO'^'G 
at f = 19.3ms. The figure shows a volume rendering of the mag- 
netic field strength (blue - green), the value of the angular veloc- 
ity, Q., at three slices parallel to the coordinate axes (red - yel- 
low), and stream lines of the velocity field in a co-rotating frame 
(black lines). In the right half of the figure, one can identify one 
large convective cell, whereas the features of, e.g., the velocity 
field in the left half are of considerably smaller scale. The mag- 
netic field is comparably weak inside the convective cell, and 
strongest in the flux tubes both at the boundaries of the convec- 
tive cell and in the left part. The colored arrows point into the 
same coordinate directions as in Fig. [19] 



4.4.3. Magnetic fields witli zero net flux 

Many results discussed above also apply analogously to models 
where the initial magnetic field has a vanishing net flux through 
the surface of the computational box. We simulated models with 
Q - 1900 s"' andoQ = -1.25, and find a basically axisymmetric 
growth of channel flows, which decay due to resistive instabili- 
ties at a level consistent with axisymmetric models. The late evo- 
lution of the models is dominated by turbulent fields. Contrary 
to uniform-field models, there is no second phase of channel ac- 
tivity. Hence, the mean magnetic energies and Maxwell stresses 
do not fluctuate violently in the saturated state, and - even when 
no velocity damping is applied - the maximum values of, e.g., 
IMct^I, are reached at MRI termination. Later on \Mj^^\ decreases 
by a factor of a few, and stays roug hly constant subsequent ly. 

Performing similar simulations lLesur & Ogilvi3 ( l2008h pro- 
posed a non-linear dynamo that balances the dissipation of the 
magnetic energy in MRI models with zero net flux. In the tur- 
bulent saturated state, they identified large-scale spatially (over 
a sizable fraction of the box) and temporally (over several rota- 
tional periods) coherent patterns of the toroidal magnetic field. 
To study this process, we looked for similar patterns in our mod- 
els. 

An example is shown in Fig. [30] for a model with b'^ - 
2 X lO'"' G simulated in a box of Ikm^ on a grid of 50^ zones 
applying no velocity damping. The figure shows the mean value 
of the toroidal field component bf^ (i.e. b'^ averaged over m and 
(f>) as a function of z and time. For t < 14ms the early channels 
flows can be identified in which the magnetic field grows. At 
t * 14ms, the channels are disrupted, and the growth of the field 
seizes. In the saturated state that follows, the mean (vertical) size 
of the structures is larger: at any time, we find only a few (typi- 
cally two) regions of opposite field polarity (blue and red), which 




-0.2 



-0.4 



Fig. 30. Evolution of large-scale coherent patterns exhibited by 
the mean (i.e., averaged over a plane at constant z-coordinate) 
toroidal field component of a 3D model having an initial mag- 
netic field of zero flux and strength b'^^ - 2x lO'^^G. The model 
was simulated in a box of Ikm^ covered by 50^ zones with 
no veloci ty damping applied. This figure is similar to Fig. 1 of 
iLesur & Qgilvie (2008.) . 



remain stable for a few rotational periods (Q~' =a . 53ms) . Thus, 
we make similar observations as iLesur & Ogilvid ( l2008h do for 
their models. 

We also simulated a few 3D models with zero-flux fields 
in the mixed regime. The results are analogous to those ob- 
tained if the initial fields are uniform(see, e.g., Fig.[3Ti). The 
MRI growth rates are similar to those of the corresponding 2D 
models, and the mean Maxwell stress in the saturated state is 
somewhat smaller. In the saturated state, both the entropy (up- 
per panel) and the angular velocity profiles become flat rapidly. 
The spatially and temporally coherent large-scale structures in 
the magnetic field are even more pronounced than in MSI mod- 
els (compare the lower panel of Fig.[3T]with Fig. [30]). They con- 
sist of two large regions characterized by an opposite sign of 

which persist for the entire simulation of the saturated state 
(a: 140ms), subject only to a slow drift in vertical direction. The 
implications of this behavior for the presence and properties of 
a non linear dynamo of the type proposed by .Lesur & OgilviS 
(I2OO8I) remain to be explored. 



5. Summary and conclusions 

We have studied the possible amplification of seed perturbations 
in supernova cores by the magneto-rotational instability. If the 
MRI grows on dynamically relevant time scales (a few tens of 
milliseconds), it can lead to MHD turbulence and efficient trans- 
port of angular momentum. Because the growth of the magnetic 
field and the associated Maxwell stresses is exponential in time, 
the MRI is one of the most promising mechanisms to amplify the 
- most likely weak - magnetic field of the supernova progenitor 

up to dynamically rel evant strengths. 

As pointed out bv lAkivama et al.l (l2003l) . the conditions for 
the instability are fulfilled in typical post-collapse supernova 
cores. Under the assumption that the MRI converts most of 
the energy contained in differential rotation into magnetic en- 
ergy, these authors predicted saturation fields of approximately 
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Fig. 31. Evolution of a mixed-regime 3D model with zero net 
flux {bl = 2x lO^^G), an entropy gradient d^S = -0.038km 
and a rotation law given by f2 = 1900 s"' and uq - -1.25. 
The upper panel shows the entropy of the model averaged over 
planes of constant tir as a function of radius and time. The lower 
panel shows, similarly to Fig.[30l the average of b"^ over planes 
of constant z as a function of z and time. 



lO'^G. This prediction derived from a semi-analytic analysis 
and ID simulations can only be confirmed by detailed multi- 
dimensional numerical simulations. The reliability of global 
simulations of the entire core, however, is limited due to the ne- 
cessity to resolve accurately small length scales (a few meters, 
at most) leading to impracticable computational costs. 

Traditionally, the MRI is studied in great detail in accretion 
discs, i.e., in systems dominated by Keplerian rotation. Because 
typical post-collapse supernova cores differ from these systems 
in many respects, e.g., by the importance of the thermal strati- 
fication and the sub-Keplerian rotation, we investigated the in- 
stability under more general physica l condi ti ons. Analyzi ng the 
MRI dispersion relation of iBalbusI (1 19951) : lUrpinI (IT996) . we 
identified the regimes of the instability relevant to supernova 
cores. 

We distinguish between Alfven and buoyant modes of the 
MRI. The former ones are generalizations of the standard MRI 
modes, and the latter ones resemble standard convective modes. 



Buoyant modes are unstable only in systems dominated by a 
negative entropy gradient, whereas Alfven modes prevail if dif- 
ferential rotation is the main agent of the instability. Whereas 
Alfven modes are rapidly amplified only for a small range of 
wave numbers, buoyant modes grow at essentially the same rate 
for a wide range of wave numbers, k < A;max- 

We have identified six regimes of the MRI depending on the 
ratio of the entropy and angular velocity gradient. These MRI 
regimes and their properties can be summarized as follows: 

1 . Sufficiently large positive gradients of the angular velocity or 
of the entropy define the stable regime with oscillatory rather 
than growing modes. 

2. For sufficiently strong differential rotation and small entropy 
gradients (or small buoyancy frequencies), we find the shear 
regime, corresponding to the hydrodynamic shear instability. 

3. If negative entropy gradients dominate the system, it is lo- 
cated in the convective regime, which resembles ordinary 
hydrodynamic convection potentially modified in the non- 
linear phase by the presence of a magnetic field. 

4. A small degree of differential rotation (e.g., Keplerian) and a 
small entropy gradient (if present at all) are the conditions for 
the magneto-shear (MSI) regime, well studied for accretion 
discs. 

5. When fast (nearly) rigid rotation suppresses convection, its 
stabilizing effect can be overridden by a weak magnetic field, 
giving rise to magneto-buoyant (MBI) modes. This regime is 
only encountered in axisymmetric flows as rotation can sta- 
bilize only axisymmetric modes of convection, i.e., in three 
dimensions convection may grow faster than the MBI. 

6. Finally, a mixed regime exists which shares properties of all 
unstable regimes listed above. 

To substantiate our stability analysis, we performed a set 
of more than 200 models of semi-global high-resolution simu- 
lations of the MRI in simplified models of post-bounce cores. 
Our novel semi-global simulations combine elements of both 
global and local simulations by taking into account the pres- 
ence of global background gradients and by providing high local 
spatial resolution. In particular, we employed t he shearing-dis c 
boundary conditions proposed by Klahr & Bod enheimeii (l2003h . 
which allow for the treatment of global gradients of, e.g., den- 
sity or entropy, and studied the influence of a thermal stratifica- 
tion on the MRI assuming various (radial) entropy profiles. The 
presence of gradients constitutes an important dilTerence of our 
setting from that of accretion discs. We used a newly developed 
Eulerian high-resolution MHD code to evolve the flow in a com- 
putational box having an edge length of a few kilometers. The 
box was located in the equatorial plane of the core at a distance 
of 15 km. The initial data were computed assuming hydrostatic 
equilibrium of differentially rotating matter described by a sim- 
plified equation of state. The gas in the box was endowed with a 
weak initial magnetic field of different topology and strength. In 
most of the simulations, the magnetic field of the progenitor had 
a strength of approximately 10'"G. We neglected the effects of 
neutrino radiation and assumed an ideal MHD flow. 

Computed under the assumption of rotational equilibrium, 
our background models differ from realistic supernova cores 
in many respects. One of the most important shortcomings is 
the neglect of the large-scale accretion flow through the post- 
shock region onto the proto-neutron star. From a physical point 
of view, such flows, whose influence on the MRI has not been 
studied previously, add to the already considerable complexity 
of the problem considered here, possibly allowing for additional 
regimes and modifying growth and saturation of the instability. 
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Technically, the proper treatment of the boundary conditions re- 
quired for a correct modeling of such flows in non-global simu- 
lations is rather involved. Thus, we have decided to postpone the 
study of the MRI in the presence of a large-scale accretion flow 
and focus on systems without overall radial motions. 

The main results of our simulations are agree well with both 
our mode analysis and with local simulations of the MRI in ac- 
cretion discs. They also confirm the estimates of lAkivama et alj 
(l2003h . and they are consistent w ith the results of global MHD 
simulations of core collapse (e.g.. lObergauhnger et al.ll2006blah . 
We summarize our results as follows: 

1. Under the physical conditions considered in this study, i.e., 
in a background rotating in hydrostatic equilibrium, the MRI 
can act in supernova cores amplifying an initial magnetic 
field strongly. The growth times are approximately equal to 
the rotational period of the core, which for rapidly rotating 
cores is sufficiently fast to influence the dynamics. 

2. Due to our relatively fine numerical grids, we were able 
to resolve the fastest growing MRI modes for initial field 
strengths higher than a few lO'^G. This threshold is con- 
sider ably lower than the one of p revious global simulations 
(e.g.. lObergaulinger et al.ll2006bllal) . enabling us to probe the 
MRI in a parameter range inaccessible to global simulations. 

3. The growth of the instability is accompanied by the devel- 
opment of channel flows, predominantly radial flows of al- 
ternating direction stacked up in the z-direction. This flow 
pattern is characteristic of both axisymmetric and three- 
dimensional simulations. The width of the channels is set 
by the wave length of the fastest growing MRI modes. 

4. At MRI termination (i.e., at the end of the first exponential 
growth phase of the instability) the channels dissolve into 
a turbulent flow having a complex magnetic field topology. 
During the subsequent evolution secondary generations of 
channel flows can (re-) appear which characterize secondary 
phases of exponential growth. 

5. We identified the mechanism responsible for the breakup of 
the channel flows and MRI termination in our simulations. 
Despite the absence of a physical resistivity in our model 
equations, we find that resistive MRI instabilities of the 
tearing-mode type develop due to the finite numerical resis- 
tivity of our MHD code. A main characteristic of the channel 
flows is the presence of prominent current sheets immersed 
between layers of opposite magnetic polarity, which are un- 
stable against current-driven instabilities. Using a simplified 
model of this kind of flows, we investigated the growth rates, 
(Tr, of the resistive instabilities, and derived an approximate 
law for the scaling of cr^ with the magnetic field strength 
present in the channels, the channel width, and the grid res- 
olution. Comparing these growth rates with the MRI ones, 
we find that the MRI ceases to grow once the tearing modes 
grow faster than the MRI. Using this criterion, we are able to 
explain the dependence of the conditions (i.e., field strength. 
Maxwell stresses) at MRI termination on, e.g., the initial 
field strength, the grid resolution, and the initial rotation pro- 
file. 

Strictly speaking, there should be no reconnection without 
physical resistivity, and the behavior of a magnetized ideal 
fluid subject to numerical resistivity may be quite differ- 
ent from that of a fluid having a large but finite conductiv- 
ity. Consequently, our results and their implications cannot 
replace a rigorous treatment of MRI growth in supernova 
cores with a non-ideal MHD model. In particular, we have 
to be careful when drawing conclusions for the MRI in non- 



ideal plasmas. Nevertheless, our results provide some qual- 
itative insight into the basic processes of MRI saturation, 
highlighting the importance of tearing-mode-like instabili- 
ties. Quantitative conclusions, as e.g., the scaling laws for the 
field strengths and the Maxwell stresses at MRI termination 
as a function of the initial magnetic field strength, should 
be taken with a grain of salt. These depend on the dissipa- 
tive properties of the numerical scheme employed, which are 
likely to change when physical resistivity is considerd. 

6. The saturation phase of the MRI differs considerably be- 
tween axisymmetric and unrestricted 3D models, and be- 
tween models having a different initial field configuration. 
In 2D the flow does not break down into small-scale turbu- 
lence, instead the channel flows merge until they form one 
pair of large-scale coherent up- and down-flows. When the 
strength of the magnetic field exceeds lO'^G, the rotational 
profile is modified within a few tens of milliseconds. 

7. Axisymmetric models having an initial magnetic field with 
a vanishing net flux through the computational box become 
turbulent after a growth phase dominated by channel flows. 
The saturation fields are considerably smaller than lO'^G. 

8. The previous finding also holds for 3D models. Turbulence 
develops, but a spontaneous reorganization of the flow may 
lead to a re-appearance of channel modes, resulting in 
Maxwell stresses comparable to those found for axisymmet- 
ric models. In models which do not develop late-stage chan- 
nel flows, field strengths up to several 10 '"^G are encoun- 
tered. The field is predominantely toroidal. The extent of 
the late-time channel activity depends on the development 
of secondary (parasitic) instabilities, both flow-driven (e.g., 
Kelvin-Helmholtz) and current-driven (e.g., tearing modes), 
which feed off the channel flows. The presence of these insta- 
bilities is determined to a large degree by the the aspect ratio 
of the computational box, i.e., we observe a strong depen- 
dence of the saturated state on the aspect ratio. For magneto- 
rotational core collapse, our results suggest that secondary 
instabilities are fairly efficient in suppressing coherent chan- 
nel flows during saturation. 

9. For models having an initial entropy gradient, we find an im- 
portant influence of convective stabilization or destabiliza- 
tion on the evolution of the MRI. We confirm the instabil- 
ity regimes predicted by our linear analysis with numerical 
simulations, the numerical growth rates being in accordance 
with the theoretical ones. The MRI is suppressed in con- 
vectively stable regions, the growth rates are reduced, and 
the geometry of the flow changes favoring radially less ex- 
tended patterns. In the mixed regime, convectively unstable 
regions with comparably large entropy gradients are domi- 
nated by flows similar to volume-filling hydrodynamic con- 
vection. The magnetic field is expelled from convective cells 
and accumulates near the box boundaries. We note that the 
entropy gradients required for these effects are fairly shallow 
~ 0.1km '. We confirm the existence of the MBI regime for 
axisymmetric models, whereas the same models, computed 
in 3D, experience the growth of non-axisymmetric modes. 

10. In 3D models having a zero net magnetic flux we observe 
the developme nt of large-scale coheren t field patterns similar 
those seen by Lesur & Q gilvie (2008"), despite the turbulent 
nature of the velocity and magnetic fields. We also find dif- 
ferences in the flow patterns between MSI and mixed regime 
models. The tentative connection to a non-linear dynamo op- 
erating in the models remains to investigated further. 
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These results allow us to draw a few conclusions. Firstly, 
the MRI has the potential to play an important role for the dy- 
namics of supernova explosions, at least for relatively fast rotat- 
ing progenitors. The details of the evolution of the MRI depend 
crucially on the properties of the core, in particular its thermal 
stratification. This makes the study of the MRI in supernovae a 
subject of its own, related to the MRI in accretion discs but also 
quite different from it. Hence, there is a need for more inves- 
tigations focusing on MRI properties specific to core collapse 
supernovae. 

While local (or semi-global) simulations can yield interest- 
ing results regarding the physics of the MRI, several important 
aspects can only be addressed by global modeling. The detailed 
dependence of the geometry of the magnetic field at saturation, 
e.g., may depend strongly on the global dynamics and on the po- 
sition of the box inside the core. Our simulations did not account 
for any of these factors: the background was in hydrostatic equi- 
librium, and we simulated models only in the equatorial region. 
Since the field geometry is of crucial importance for the global 
dynamics, e.g., for the generation of jet-like outflows in collap- 
sars, conclusions on the dynamic influence of the MRI based on 
local simulations cannot be drawn easily. They require global 
models. 

Additional investigations are also required to study the inter- 
play of large-scale flows, e.g., the accretion of matter onto the 
proto-neutron star and the influence of neutrino radiation. Apart 
from modifying the regimes of the instability discussed here, the 
former process may allow for new instability regimes; in particu- 
lar the combination of the MRI and the standing accretion shock 
instability may have interesting consequences for the dynamics 
of the explosion. While neutrino heating and cooling may, due 
to a net transport of entropy, lead to a stratification qualitatively 
similar to that of some models in this study with non-vanishing 
entropy gradients, it is too early to speculate on how neutrino 
radiation affects the MRI. 

The inclusion of the MRI and its effects in global simula- 
tions requires a considerably careful treatment. The currently 
used approach of artificially enhancing the initial field strength 
by a constant factor is questionable. On the other hand, finding 
a better prescription relies on unraveling the dependence of sat- 
urated MRI driven turbulence on the different parameters of the 
system, as e.g., the rotation law, the thermodynamic conditions, 
and probably also the neutrino transport. 

With our current simulations we are unfortunately not yet 
able to go beyond the stage of a qualitative proof of principle and 
to address important open questions of the MRI in core collapse 
supernovae. This would require additional 3D high-resolution 
non-ideal MHD simulations covering a large parameter space of 
possible rotation profiles, thermal stratifications, and magnetic 
field geometries. We are planning to address these issues in fu- 
ture work. 
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Appendix A: Growth rates of resistive instabilities 

Our simulations of the MRI indicate that the termination of the 
growth of the instability is determined, at least partially, by resis- 
tive instabilities of the tearing-mode type. Although there exist 
deta iled inves tigations of this kind of resistive instabilities (see, 
e.g.. lBiskampli2000,) . the application of the results to our study 
is hampered by the completely different type of dissipative ef- 
fects we are facing here: all previous results hold for instabilities 
due to physical resistivity, whereas our ideal MHD simulations 
are affected by numerical resistivity, only. Hence, we had to de- 
termine the growth rates of resistive instabilities from numerical 
experiments without referring to analytic results - although, as 
we will see, there exist certain similarities. 

We simulated the evolution of two-dimensional current-sheet 
models on a Cartesian grid [jcq, x\\ x [yo.yi] with x zones 
imposing periodic boundary conditions. The fluid is described 
by an ideal-gas equation of state with an adiabatic index F = 4/3. 
We used initial data mimicking MRI-generated channel flows: 
the initial magnetic field varies sinusoidally in a gas of constant 
density po and pressure Pq, 



b - /jgsm , 

and having a velocity in x-direction given by 



V — Vr, sm ■ 



Iniy - yo) - 7r/2 



(A.l) 
(A.2) 

(A.3) 



Here, a denotes the initial width of the flux sheet, and is equal 
to one half of the Alfven velocity corresponding to b-^. The 
presence of the (shear-free!) initial x- velocity is not essential, as 
it changes the growth rates of the instability only little. However, 
as we observe this kind of a velocity in channel flows, we have 
included it in our simulations. 

We perturbed the sheet by a smaU random j- velocity (10^^ x 
c^.). The parameters of the model are chosen to mimic the situa- 



tions encountered in MRI simulations (see Tab. lA.lb . To isolate 
the dependence of the growth rates on the physical and numeri- 
cal parameters, we varied the initial magnetic field strength, b^, 
the initial density, po, the width of the current sheet, a, and the 
grid resolution, 6x - (xi - xo)/m^. We chose the grid resolution 
such that the current sheet is covered by 12 to 100 zones. The 
initial pressure is Po = xp^. 

We show one typical result for the evolution of our models 
in Fig. lA.ll After a short initial phase, the transverse magnetic 
energy density (green line) grows roughly exponentially as tear- 
ing modes develop. Initially the growth rate is approximately 
constant, but it increases by a factor of ~ 4 towards saturation. 
Simultaneously, the jc-component of the magnetic field decreases 
strongly until it is of similar strength as the y-component. At this 
point, the coherent current sheets are completely disrupted by 
the resistive instability. 

We determined estimates of the growth rates of the resistive 
instability using the time derivative of the transverse magnetic 
energy density, i.e., the time derivative of the magnetic energy 
density corresponding to the y-component of the magnetic field. 
To find a scaling relation of the form 



^ (4) ('^s)''^ (SxY^ (aY" 



we define the following functions of the growth rate cr. 



fcA 
fa 



cr{ciy''\csrHSxr^' ■ 
cr(csr^Har^'(6x)-^\ 



(A.4) 

(A.5) 
(A.6) 
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Table A.l. Parameters of the 2D current sheets simulated to de- 
termine the growth rates of resistive instabilities. The columns 
give from left to right the edge length of the square simulation 
box, the number of grid zones per dimension, the initial density 
(in units of lO''' g/cm"^), the sound speed (in units of 10^ cm/s), 
the magnetic field (in units of 10'^ G), the wavelength of the ini- 
tial field, and the growth rate, respectively. 



L 




P0;13 


CS;9 




a 


a 


[km] 




[g cm^'l 
L= 1 


cm s"' 
1 1 


[G] 


[m] 


[ms-'l 


0.25 


50 


2.5 


3.1 


8 


125 


0.45 


0.25 


50 


2.5 


3.1 


16 


125 


1.6 


0.25 


50 


2.5 


3.1 


32 


125 


6.0 


0.25 


50 


2.5 


3.1 


64 


125 


32 


0.25 


50 


2.5 


3.1 


16 


62.5 


10 


0.25 


50 


2.5 


3.1 


32 


62.5 


23 




1 fin 


9 S 

Z.J 


^ 1 


1 ft 


1 95 

IZJ 


n 6n 

u.uu 




1 nn 

lUU 


9 S 

Z.J 


^ 1 


^9 

JZ 


1 95 

IZJ 


J. i 


0.25 


100 


2.5 


3.1 


64 


125 


12 




100 


2.5 


3.1 


128 


125 


43 


U.ZJ 


1 no 


Z.J 


1 1 

J. I 


iO 


A9 s 






1 on 

iUU 


Z.. J 




'?9 

JZ 














64 
ut 




SO 


0.25 


100 


2.5 


3.1 


128 


62.5 


110 




100 


0.025 


1.4 


^ 9 
J .z 


UZ. J 


9S 

Z J 




1 nn 

IUU 


n n9'^ 

U.UZ J 


D. L 


^ 9 
J .z 


69 5 

UZ. J 


i J 


n 9s 


1 nn 


n n9's 

U.UZ. J 


U. i 


^ 9 
J .z 


69 5 

UZ. J 


8 n 

o.u 




1 nn 

iUU 


9 S 

Z.J 


1 d 


^9 

JZ 


69 5 

UZ. J 


Z J 




1 nn 

iUU 


9 S 

Z.J 


^ 1 
J. 1 


^9 

JZ 


69 5 

UZ. J 


it 


0.25 


100 


2.5 


1.4 


6.1 


62.5 


7.8 


f) 95 


inn 

IUU 


95n 

ZJU 


u.u 


^9n 

JZU 


69 5 

UZ. J 


7.4 


0.25 


100 


250 


3.1 


320 


62.5 


15 


0.25 


100 


9 'S 

Z.J 


12.4 


64 

Ut- 


69 S 

UZ. J 


1 6 
i u 


0.25 


100 


2.5 


3.1 


4 


31.25 


1.2 


0.25 


100 


2.5 


3.1 


g 


31.25 


5.9 


f) 9S 


inn 

IUU 


2.5 


3.1 


16 


31.25 


20 


0.25 


100 


2.5 


3.1 


32 


31.25 


53 


0.25 


100 


2.5 


3.1 


64 


31.25 


110 


0.25 


200 


2.5 


3.1 


32 


125 


0.9 


0.25 


200 


2.5 


3.1 


64 


125 


4.6 


0.25 


200 


2.5 


3.1 


128 


125 


19 


0.25 


200 


2.5 


3.1 


32 


62.5 


5.4 


0.25 


200 


2.5 


3.1 


64 


62.5 


24 


0.25 


200 


2.5 


3.1 


128 


62.5 


76 


0.25 


200 


2.5 


3.1 


16 


31.25 


7.9 


0.25 


200 


2.5 


3.1 


32 


31.25 


28 


0.25 


200 


2.5 


3.1 


64 


31.25 


95 


0.25 


200 


2.5 


3.1 


128 


31.25 


250 


0.25 


200 


2.5 


3.1 


4 


15.625 


2.5 


0.25 


200 


2.5 


3.1 


16 


15.625 


40 


0.25 


200 


2.5 


3.1 


64 


15.625 


260 


0.5 


100 


2.5 


3.1 


32 


250 


1.2 


0.5 


100 


2.5 


3.1 


64 


250 


5.6 


0.5 


100 


2.5 


3.1 


128 


250 


20 


0.5 


100 


2.5 


3.1 


16 


125 


1.9 


0.5 


100 


2.5 


3.1 


32 


125 


7.5 


0.5 


100 


2.5 


3.1 


64 


125 


26 


0.5 


100 


2.5 


3.1 


8 


62.5 


3.0 


0.5 


100 


2.5 


3.1 


16 


62.5 


10 


0.5 


100 


2.5 


3.1 


32 


62.5 


24 


0.5 


100 


2.5 


3.1 


64 


62.5 


61 



and adjust the exponents jp^, ys. Ta^ and to determine the scal- 
ing of cr with the respective parameters. Our preferred set of 
scaling exponents is 

TA = 1.75, (A.7) 



^ni : : 
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Fig. A.l. Evolution of the magnetic energy density of a current- 
sheet model simulated on a square computational grid (edge 
length 0.25 km; 200 x 200 zones). The model has an initial 
density po — 2.5 x lO'^gcm"^, and an initial field strength 
/7y = 1.6 X lO'^G. The wavelength of the initially sinusoidally 
varying magnetic field is a = 3 1 .25m. The black solid line and 
the green dashed line show the magnetic energy density corre- 
sponding to the X and the y component of the magnetic field, 
respectively. 

rs = -0.75, (A.8) 
To = -2, (A.9) 
75 = 1, (A. 10) 

which implies the following scaling law: 




(A.ll) 



We demonstrate the quality of the fit parameters in Fig. lA.2l 
showing fcA as a function of the initial Alfven velocity (left 
panel), and fa as a function of the width of the current sheet 
(middle panel) and of the grid resolution (right panel). None of 
the groups of models representing the variation of one parameter 
(distinguished by a different color and symbol in the figure) ex- 
hibits a strong trend with ca, a or 5x, i.e., our scaling exponents 
provide an adequate fit to the data. 

Due to a relatively large scatter in the growth rates, the scal- 
ing relation, Eq. ( IA.4b . should not be taken too literally. We note, 
however, that our computed rates are compatible with those of 
the resistive instabilities (for fields of similar strength) in MRI 
models, i.e., approximately one millisecond for ~ lO'^G. An 
additional dependence of cr on the domain size, L, cannot be ex- 
cluded, but we did not examine this possibility any further. For 
small a and coarse grids our scaling formula tends to overesti- 
mate the growth rate for large initial Alfven velocities, and for 
sound speeds much larger than the Alfven velocity the rates de- 
pend more strongly on the sound speed than predicted by our 
formula. Because both situations do not apply to our MRI mod- 
els, we did not pursue these issues any further Our scaling law 
loses its validity, if the Alfven velocity exceeds the sound speed. 
Thus, we excluded two respective models in the derivation of 
our scaling relation. 

Bearing in mind the uncertainties regarding the physical 
meaning of a purely numerical resistivity and the precise val- 
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ues of the scaling we may try to interpret our result summa- 
rized in Eq. dA.llt . As the product of the Alfven velocity and 
the grid spacing, cp^5x, defines an effective resistivity, we may 
conclude that the growth time of the instability is set by the time 
scale for resistive diffusion across the width of a current sheet, 
Tr = HcaSx), modified by the ratio of the sound speed and the 
Alfven velocity. This interpretation has the nice property that it 
is consistent with the fact that the magnetic Reynolds number is 
proportional to the grid resolution. 
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Appendix B: List of models 

In this appendix we provide a list of the models: 

Tab. lB.ll : This table contains a list of 2D models having a pos- 
itive entropy gradient. Their initial rotation profile is given 
by Oo = 1900 s"' and an - -1.25, and their initial magnetic 
field is uniform. 

Tab. IB. 21 : This table contains a list of 2D models having a neg- 
ative entropy gradient. The models rotate initially rigidly or 
differentially. The initial magnetic field is uniform. 

Tab. |B.3| and Tab. |B.4| : This table contains a list of 3D models 
having different initial magnetic field strengths, entropy gra- 
dients, initial rotation laws, and simulated in computational 
boxes of various size with both types (with and without ve- 
locity damping) of radial boundary conditions. 
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Fig. A.2. Dependence of the growth rate of resistive instabilities, cr, on various parameters of 2D current sheet models. The panels 
show fcA as a function of the initial Alfven velocity (c;^; left panel), and /„ as a function of the width of the current sheet (a; 
middle panel) and of the grid resolution (6x; right panel). The left panel shows groups of models with different a: 250 m (pink 
plus sign), 125 m (blue diamond), 62.5 m (green triangle), 32.25 m (orange square), and 15.625 m (red cross). The other two panels 
show groups of models with different Alfven velocity (in units of lO^cm s"'): 0.8 (black asterisk), 1.6 (pink plus sign), 3.2 (blue 
diamond), 6.4 (green triangle), 12.8 (orange square), and 25.6 (red cross) respectively. 
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Table B.l. List of 2D models having a positive entropy gradient and an initial rotation profile given by - 1900 s ' and an = 
-1.25. The columns give (from left to right) the number of grid zones nij^xm, (the box has an edge length L^^ - L,- 1km), the type 
of boundary condition which was applied (d: velocity damping; p: periodic), the adiabatic index of the gas, F (models computed 
with an ideal-gas equation of state instead of the hybrid one are marked "id"), the initial entropy, 5o, and the radial entropy gradient, 
dj^S . The next three columns give the parameter C determining the instability regime (see Eq. (l24li). the theoretical growth rate, o-jh, 
and the strength of the initially uniform magnetic field, Z?^. The last two columns list the numerical growth rate, cr, and the value of 
the Maxwell stress component M^™ at MRI termination. 
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Table B.2. List of 2D models having a negative initial entropy gradient. The columns give (from left to right) the number of grid 
zones m^j x m. (the box has an edge length - L, = 1km), the type of boundary condition which was applied (d: velocity 
damping; p: periodic), the rotation law ("d": differential rotation with f2o = 1900 s-' and an = -1.25; "r"" rigid rotation with an 
angular velocity of Q), the adiabatic index of the gas, F, the initial entropy. So, the radial entropy gradient, d^S , and the strength of 
the initially uniform magnetic field, bfy The next columns give the three quantities H^IQ^y N^/i^l and C determining the instability 
regime (see Eq. (|24]|). followed by the theoretical growth rate, cTth, and the type of the instability. The last two columns list the 
numerical growth rate, cr, and the value of the Maxwell stress component M^™ at MRI termination. 
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Table B.3. List of 3D models. The first column (from left to right) gives the geometry of the initial field ("U": uniform; "V" zero- 
flux). The next two columns show the box size (L^ X X L^) and the number of grid zones (Wro^ X X m^). The next four columns 
Ust the rotation law ("d": differential rotation with Qq = 1900 s~^ and an = -1.25; "r'^" rigid rotation with an angular velocity of 
Q), the entropy gradient, {sq - 0.2), the strength of the initial magnetic field, and the type of boundary condition which was 
applied (d: velocity damping; p: periodic). The remaining four columns give the growth rate of the MRI, cr, and the Maxwell stress 
component M^T* at MRI termination, its maximum value, and its time averaged value, respectively. 
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Table B.4. Continuation of Tab. |B3] For the MBI model with /jg - lOG, we do not provide values for the Maxwell stress, because 
the magnetic field behaves similar to a purely passive field. 
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